Find significantly different loci/regions between years using shuffling to create null distributions

source("../Rscripts/BaseScripts.R")
require(data.table)
require(plyr)
require(RColorBrewer)

1 Theta null distribution

1.1 Create persite theta.pest.gz files in angsd

```bash

/home/jamcgirr/apps/angsd/misc/thetaStat  print /home/ktist/ph/data/angsd/theta/PWS91_maf00.thetas.idx | gzip > /home/ktist/ph/data/angsd/theta/PWS91_maf00_thetas.pestPG.gz

# Run printTheta.sh at farm

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->



## Run R scripts at farm to create null distribution of genome-wide thetas

Run angsd_theta_siteshuffle_null.sh at farm, which runs Pi_shuffle_pws.R 
 - Takes long time, so create a script for each pop.year combination and run separately 

Output from theta shuffling results are in Data/shuffle/theta.siteshuffle.50000.PWS91_PWS96.csv.gz


## Calculate p-values for differences in theta, pi, and Tajima's D using the results from shuffling 


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjpbIiNmcm9tIGNvZEV2b2wgYW5nc2RfdGhldGFfc2l0ZXNodWZmbGVfbnVsbF9zdGF0cy5SIiwiIiwiIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIiwiIyBsb2FkIGZ1bmN0aW9ucyIsIiMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyIsInJlcXVpcmUoZGF0YS50YWJsZSkiLCJyZXF1aXJlKHBseXIpIiwicmVxdWlyZShnZ3Bsb3QyKSIsInJlcXVpcmUoUkNvbG9yQnJld2VyKSIsIiIsIiIsIiMjIyMjIyMjIyMjIyMjIyMjIyMjIyIsIiMgcmVhZCBpbiBhbmQgcHJlcCBkYXRhIiwiIyMjIyMjIyMjIyMjIyMjIyMjIyMjIiwiIiwieWVhcnM8LWMoXCJQV1M5MVwiLFwiUFdTOTZcIixcIlBXUzA3XCIsXCJQV1MxN1wiKSIsImNvbWI8LXQoY29tYm4oeWVhcnMsIDIpKSIsIiIsImNocm1heCA8LSBmcmVhZCgnLi4vRGF0YS9uZXdfdmNmL2Nocl9zaXplcy5iZWQnKSIsImNocm1heDwtY2hybWF4WywtMl0iLCJjb2xuYW1lcyhjaHJtYXgpPC1jKFwiY2hyXCIsIFwibGVuXCIpIiwiY2hybWF4JHN0YXJ0PC1jKDAsY3Vtc3VtKGNocm1heCRsZW4pWzE6KG5yb3coY2hybWF4KS0xKV0pIiwiIiwiY2hybWF4JGVuZDwtY3Vtc3VtKGNocm1heCRsZW4pIiwiY2hybWF4JG1pZGRsZTwtKGNocm1heCRlbmQtY2hybWF4JHN0YXJ0KS8yK2Nocm1heCRzdGFydCIsIiIsIiNzZXRrZXkoY2hybWF4LCBjaHIpIiwiIiwiI0Z1bmN0aW9ucyB0byBjYWxjdWxhdGUgcC12YWx1ZXMgZnJvbSBjb2RFdm9sIiwiY2FsY3BHIDwtIGZ1bmN0aW9uKHRoZXRhY2hhbmdlLCBudWxsKXsgIyBmb3IgaW5jcmVhc2VzIGluIHRoZXRhIiwiICByZXR1cm4oKHN1bShudWxsID4gdGhldGFjaGFuZ2UpKzEpLyhsZW5ndGgobnVsbCkrMSkpICMgZXF1YXRpb24gZnJvbSBOb3J0aCBldCBhbC4gMjAwMiBBbSBKIEh1bSBHZW4iLCJ9IiwiY2FsY3BMIDwtIGZ1bmN0aW9uKHRoZXRhY2hhbmdlLCBudWxsKXsgIyBmb3IgZGVjcmVhc2VzIGluIHRoZXRhIiwiICByZXR1cm4oKHN1bShudWxsIDwgdGhldGFjaGFuZ2UpKzEpLyhsZW5ndGgobnVsbCkrMSkpICMgZXF1YXRpb24gZnJvbSBOb3J0aCBldCBhbC4gMjAwMiBBbSBKIEh1bSBHZW4iLCJ9IiwiIiwiIiwiY2hfY29scyA8LSBicmV3ZXIucGFsKDQsICdQYWlyZWQnKVtyZXAoMToyLDEzKV0iLCJEYXRhbGw8LWRhdGEudGFibGUoKSIsImZvciAocCBpbiAxOiBucm93KGNvbWIpKXsiLCIgICAgIyBtYXggdGhldGEgcGVyIGdlbm9tZSBmcm9tIHJlc2h1ZmZsaW5nIChhbGwgc2l0ZXMpIGZyb20gYW5nc2RfdGhldGFfc2l0ZXNodWZmbGVfbnVsbC5yIiwiICAgIG51bGw8LWZyZWFkKHBhc3RlMCgnLi4vRGF0YS9zaHVmZmxlL3RoZXRhLnNpdGVzaHVmZmxlLjUwMDAwLicsIGNvbWJbcCwxXSxcIl9cIixjb21iW3AsMl0sJy5jc3YuZ3onKSkiLCIgICAgIiwiICAgICN1cHBlciBhbmQgbG93ZXIgOTUlIiwiICAgIG51bGxbLCAuKHRXZF9sOTUgPSBxdWFudGlsZShtaW50V2QsIDAuMDUpLCB0V2RfdTk1ID0gcXVhbnRpbGUobWF4dFdkLCBwcm9icyA9IDAuOTUpLCIsIiAgICAgICAgICAgICAgICB0UGRfbDk1ID0gcXVhbnRpbGUobWludFBkLCAwLjA1KSwgdFBkX3U5NSA9IHF1YW50aWxlKG1heHRQZCwgcHJvYnMgPSAwLjk1KSwiLCIgICAgICAgICAgICAgICAgdERkX2w5NSA9IHF1YW50aWxlKG1pbnREZCwgMC4wNSksIHREZF91OTUgPSBxdWFudGlsZShtYXh0RGQsIHByb2JzID0gMC45NSkpXSIsIiIsIiAgICAjYXNzaWduKHBhc3RlMChcIm51bGwuXCIsY29tYltwLDFdLFwiX1wiLGNvbWJbcCwyXSksIG51bGwpICAgICIsIiAgICAiLCIgICAgIyBzbGlkaW5nIHdpbmRvd3MgdGhldGEgY2hhbmdlIChHQVRLIHNpdGVzKSBmcm9tIGFuZ3NkX3RoZXRhX3NpdGVzaHVmZmxlX251bGwuciIsIiAgICBkYXQ8LWZyZWFkKHBhc3RlMCgnLi4vRGF0YS9zaHVmZmxlL3RoZXRhX2NoYW5nZV9yZWdpb25fNTAwMDAuJywgY29tYltwLDFdLFwiX1wiLGNvbWJbcCwyXSwnLmNzdi5neicpLCBkcm9wID0gMSkiLCIgICAgZGF0Wyxwb3A6PXBhc3RlMChjb21iW3AsMV0sXCJfXCIsY29tYltwLDJdKV0iLCIgICAgIiwiICAgIGRhdDwtbWVyZ2UoZGF0LCBjaHJtYXhbLGMoXCJjaHJcIixcInN0YXJ0XCIpXSwgYnkueD1cIkNocm9tb1wiLCBieS55ID0gXCJjaHJcIikiLCIgICAgZGF0WywgUE9TZ2VuIDo9IFdpbkNlbnRlciArIHN0YXJ0XSIsIiAgICBkYXRbLHN0YXJ0IDo9IE5VTExdICNyZW1vdmUgc3RhcnQiLCIgICAgIiwiICAgICNjYWxjdWxhdGUgcC12YWx1ZXMiLCIgICAgIzEuIHRoZXRhVyBsb2NpIiwiICAgIGRhdFt0V2QgPiAwLCB0V2QucCA6PSBjYWxjcEcodFdkLCBudWxsJG1heHRXZCksIGJ5ID0gLihDaHJvbW8sIFdpbkNlbnRlcildICMgdGhldGFXICBsb2NpIiwiICAgIGRhdFt0V2QgPD0gMCwgdFdkLnAgOj0gY2FsY3BMKHRXZCwgbnVsbCRtaW50V2QpLCBieSA9IC4oQ2hyb21vLCBXaW5DZW50ZXIpXSIsIiAgICAjMi4gdGhldGEgcGkiLCIgICAgZGF0W3RQZCA+IDAsIHRQZC5wIDo9IGNhbGNwRyh0UGQsIG51bGwkbWF4dFBkKSwgYnkgPSAuKENocm9tbywgV2luQ2VudGVyKV0gIyB0aGV0YSBwaSIsIiAgICBkYXRbdFBkIDw9IDAsIHRQZC5wIDo9IGNhbGNwTCh0UGQsIG51bGwkbWludFBkKSwgYnkgPSAuKENocm9tbywgV2luQ2VudGVyKV0iLCIgICAgIiwiICAgICNUYWppbWEncyBEIiwiICAgIGRhdFt0RGQgPiAwLCB0RGQucCA6PSBjYWxjcEcodERkLCBudWxsJG1heHREZCksIGJ5ID0gLihDaHJvbW8sIFdpbkNlbnRlcildICMgdGFqaW1hJ3MgRCIsIiAgICBkYXRbdERkIDw9IDAsIHREZC5wIDo9IGNhbGNwTCh0RGQsIG51bGwkbWludERkKSwgYnkgPSAuKENocm9tbywgV2luQ2VudGVyKV0iLCIiLCIgICAgd3JpdGUuY3N2KGRhdCwgZmlsZT1nemZpbGUocGFzdGUwKCcuLi9PdXRwdXQvUGkvU2h1ZmZsZS90aGV0YV9zaXRlc2h1ZmZsZV8nLCBjb21iW3AsMV0sXCJfXCIsY29tYltwLDJdLCcuY3N2Lmd6JykpKSIsIiAgICAiLCIgICAgRGF0YWxsPC1yYmluZChEYXRhbGwsIGRhdCkiLCJ9IiwiIiwid3JpdGUuY3N2KERhdGFsbCwgZmlsZT1nemZpbGUocGFzdGUwKCcuLi9PdXRwdXQvUGkvU2h1ZmZsZS90aGV0YV9zaXRlc2h1ZmZsZV9QV1Nfc3VtbWFyeS5jc3YuZ3onKSkpIl19 -->

```r
#from codEvol angsd_theta_siteshuffle_null_stats.R

###########################
# load functions
###########################
require(data.table)
require(plyr)
require(ggplot2)
require(RColorBrewer)


#####################
# read in and prep data
#####################

years<-c("PWS91","PWS96","PWS07","PWS17")
comb<-t(combn(years, 2))

chrmax <- fread('../Data/new_vcf/chr_sizes.bed')
chrmax<-chrmax[,-2]
colnames(chrmax)<-c("chr", "len")
chrmax$start<-c(0,cumsum(chrmax$len)[1:(nrow(chrmax)-1)])

chrmax$end<-cumsum(chrmax$len)
chrmax$middle<-(chrmax$end-chrmax$start)/2+chrmax$start

#setkey(chrmax, chr)

#Functions to calculate p-values from codEvol
calcpG <- function(thetachange, null){ # for increases in theta
  return((sum(null > thetachange)+1)/(length(null)+1)) # equation from North et al. 2002 Am J Hum Gen
}
calcpL <- function(thetachange, null){ # for decreases in theta
  return((sum(null < thetachange)+1)/(length(null)+1)) # equation from North et al. 2002 Am J Hum Gen
}


ch_cols <- brewer.pal(4, 'Paired')[rep(1:2,13)]
Datall<-data.table()
for (p in 1: nrow(comb)){
    # max theta per genome from reshuffling (all sites) from angsd_theta_siteshuffle_null.r
    null<-fread(paste0('../Data/shuffle/theta.siteshuffle.50000.', comb[p,1],"_",comb[p,2],'.csv.gz'))
    
    #upper and lower 95%
    null[, .(tWd_l95 = quantile(mintWd, 0.05), tWd_u95 = quantile(maxtWd, probs = 0.95),
                tPd_l95 = quantile(mintPd, 0.05), tPd_u95 = quantile(maxtPd, probs = 0.95),
                tDd_l95 = quantile(mintDd, 0.05), tDd_u95 = quantile(maxtDd, probs = 0.95))]

    #assign(paste0("null.",comb[p,1],"_",comb[p,2]), null)    
    
    # sliding windows theta change (GATK sites) from angsd_theta_siteshuffle_null.r
    dat<-fread(paste0('../Data/shuffle/theta_change_region_50000.', comb[p,1],"_",comb[p,2],'.csv.gz'), drop = 1)
    dat[,pop:=paste0(comb[p,1],"_",comb[p,2])]
    
    dat<-merge(dat, chrmax[,c("chr","start")], by.x="Chromo", by.y = "chr")
    dat[, POSgen := WinCenter + start]
    dat[,start := NULL] #remove start
    
    #calculate p-values
    #1. thetaW loci
    dat[tWd > 0, tWd.p := calcpG(tWd, null$maxtWd), by = .(Chromo, WinCenter)] # thetaW  loci
    dat[tWd <= 0, tWd.p := calcpL(tWd, null$mintWd), by = .(Chromo, WinCenter)]
    #2. theta pi
    dat[tPd > 0, tPd.p := calcpG(tPd, null$maxtPd), by = .(Chromo, WinCenter)] # theta pi
    dat[tPd <= 0, tPd.p := calcpL(tPd, null$mintPd), by = .(Chromo, WinCenter)]
    
    #Tajima's D
    dat[tDd > 0, tDd.p := calcpG(tDd, null$maxtDd), by = .(Chromo, WinCenter)] # tajima's D
    dat[tDd <= 0, tDd.p := calcpL(tDd, null$mintDd), by = .(Chromo, WinCenter)]

    write.csv(dat, file=gzfile(paste0('../Output/Pi/Shuffle/theta_siteshuffle_', comb[p,1],"_",comb[p,2],'.csv.gz')))
    
    Datall<-rbind(Datall, dat)
}

write.csv(Datall, file=gzfile(paste0('../Output/Pi/Shuffle/theta_siteshuffle_PWS_summary.csv.gz')))

1.2 Plot the results

1.2.1 Changes in Pi/Theta/D between years

## Plot the results ##
Datall<-fread('../Output/Pi/Shuffle/theta_siteshuffle_PWS_summary.csv.gz')

winsz = 5e4 

#Changes in Pi between years
Datall$Chromo<-factor(Datall$Chromo, levels=c(paste0("chr",1:26)))
Datall$pop<-factor(Datall$pop, levels=c("PWS91_PWS96","PWS91_PWS07","PWS91_PWS17","PWS96_PWS07","PWS96_PWS17","PWS07_PWS17"))
ch_cols <- brewer.pal(4, 'Paired')[rep(1:2,13)]
ggplot(Datall, aes(POSgen, tPd/winsz, color = Chromo)) + 
    geom_point(size = 0.5, alpha = 0.3) +
    facet_wrap(~pop, ncol = 1) +
    scale_color_manual(values = ch_cols, guide="none") +
    ylab('Change in pi per site')+xlab("Chromosome")+
    ggtitle("Changes in Pi")+
    scale_x_continuous(breaks=chrmax$middle, labels=1:26)+
    theme_bw()
ggsave(paste0('../Output/Pi/Shuffle/Changes_in_Pi_PWS.png'), width = 7.5, height = 9, dpi = 300)

# plot theta_Watterson change
ggplot(Datall, aes(POSgen, tWd/winsz, color = Chromo)) + 
  geom_point(size = 0.5, alpha = 0.3) +
  scale_color_manual(values = ch_cols, guide="none") +
    facet_wrap(~pop, ncol = 1) +
  ylab('Change in Wattersons theta per site')+xlab("Chromosome")+
  ggtitle("Changes in Pi")+
    scale_x_continuous(breaks=chrmax$middle, labels=1:26)+
    theme_bw()
ggsave('../Output/Pi/Shuffle/Changes_in_thetaW_PWS.png', width = 7.5, height = 9, dpi = 300)

# plot Tajima's D change
ggplot(Datall, aes(POSgen, tDd, color = Chromo)) + 
    geom_point(size = 0.5, alpha = 0.3) +
    scale_color_manual(values = ch_cols,guide="none") +
    facet_wrap(~pop, ncol = 1) +
    ylab('Change in Tajimas D per window')+xlab("Chromosome")+
    ggtitle("Changes in Tajima's D")+
    scale_x_continuous(breaks=chrmax$middle, labels=1:26)+
    theme_bw()
ggsave('../Output/Pi/Shuffle/Changes_in_TajimasD_PWS.png', width = 7.5, height = 9, dpi = 300)

1.2.2 Plot the p-values for differences in P/Theta along the genome

# plot pi p-value vs. position 
ch_cols <- brewer.pal(4, 'Paired')[rep(1:2,13)]
Datall$Chromo<-factor(Datall$Chromo, levels=c(paste0("chr",1:26)))
Datall$pop<-factor(Datall$pop, levels=c("PWS91_PWS96","PWS91_PWS07","PWS91_PWS17","PWS96_PWS07","PWS96_PWS17","PWS07_PWS17"))

ggplot(Datall, aes(POSgen, -log10(tPd.p)*sign(tPd), color = Chromo)) + 
    geom_point(size = 0.4, alpha = 0.3) +
    facet_wrap(~pop, ncol = 1) +
    scale_color_manual(values = ch_cols, guide="none") +xlab("Chromosome")+
    ylab("log10(P-value)")+
    geom_hline(yintercept = log10(0.05), linetype = 'dashed', color = 'red',size=0.3) +
    geom_hline(yintercept = -log10(0.05), linetype = 'dashed', color = 'red', size=0.3)+
    ggtitle("P-values for changes in Pi")+
    scale_x_continuous(breaks=chrmax$middle, labels=1:26)+
    theme_bw()
ggsave('../Output/Pi/Shuffle/Changes_in_Pi.siteshuffle.p-values_PWS.png', width = 7.5, height =11, units = 'in', dpi = 300)


# plot thetaW p-value vs. position (all loci)
ggplot(Datall, aes(POSgen, -log10(tWd.p)*sign(tWd), color = Chromo)) + 
  geom_point(size = 0.4, alpha = 0.3) +
  facet_wrap(~pop, ncol = 1) +
  scale_color_manual(values = ch_cols, guide="none") +xlab("Chromosome")+
    ylab("log10(P-value)")+
  geom_hline(yintercept = log10(0.05), linetype = 'dashed', color = 'red', size=0.3) +
  geom_hline(yintercept = -log10(0.05), linetype = 'dashed', color = 'red', size=0.3)+
      ggtitle("P-values for changes in Theta")+
    scale_x_continuous(breaks=chrmax$middle, labels=1:26)+
    theme_bw()
ggsave('../Output/Pi/Shuffle/Changes_in_thetaW.siteshuffle.p-values_PWS.png', width = 7.5, height =11, units = 'in', dpi = 300)

#plot Tajama's D p-value vs. position
ggplot(Datall, aes(POSgen, -log10(tDd.p)*sign(tDd), color = Chromo)) + 
    geom_point(size = 0.4, alpha = 0.3) +
    facet_wrap(~pop, ncol = 1) +
    scale_color_manual(values = ch_cols, guide="none") +xlab("Chromosome")+
    ylab("log10(P-value)")+
    geom_hline(yintercept = log10(0.05), linetype = 'dashed',  color = 'red', size=0.3) +
    geom_hline(yintercept = -log10(0.05), linetype = 'dashed', color = 'red', size=0.3)+
    ggtitle("P-values for changes in Tajima's D")+
    scale_x_continuous(breaks=chrmax$middle, labels=1:26)+
    theme_bw()
ggsave('../Output/Pi/Shuffle/Changes_in_TajimaD.siteshuffle.p-values_PWS.png', width = 7.5, height =11, units = 'in', dpi = 300)

1.3 Find the regions with signfiicant P-values (Outlier regions)

1.3.1 Pi

# Outliers

Pi_outliers<-Datall[tPd.p < 0.05,]
Theta_outliers<-Datall[tWd.p < 0.05,]
TajimaD_outliers<-Datall[tDd.p < 0.05,] #no outliers

# Summary per chromosome per population
pi<-data.frame(table(Pi_outliers$pop, Pi_outliers$Chromo))
the<-data.frame(table(Theta_outliers$pop, Theta_outliers$Chromo))
D<-data.frame(table(TajimaD_outliers$pop, TajimaD_outliers$Chromo))

# Summary per population
pi2<-data.frame(table(Pi_outliers$pop))
the2<-data.frame(table(Theta_outliers$pop))
Ds2<-data.frame(table(TajimaD_outliers$pop))


#plot PWS91-96, 96-07, and 07-17
yrs<-c("PWS91_PWS96","PWS96_PWS07","PWS07_PWS17")
col3<-brewer.pal(5,"PuRd")[c(2,3,5)]
div1<-diverging_hcl(6, palette="Blue-Red")
#reverse the color
div2<-rev(div1)

pis<-pi[pi$Var1 %in% yrs,]
pis$Var1<-factor(pis$Var1, levels=yrs)
ggplot(pis, aes(x=Var2, y=Freq, fill=Var1))+
    geom_bar(stat="identity",position=position_dodge(width=0.8))+
    scale_fill_manual(values=col3)+
    theme_classic()+
    theme(axis.text.x = element_text(angle=45, hjust=1), legend.title = element_blank())+
    ggtitle(expression(paste("Changes in ", pi)))+
    xlab('')+ylab('Number of regions with P<0.05')
ggsave("../Output/Pi/Shuffle/Pi_significant_perChrom_perPop.png", width = 8, height = 4, dpi=300) 

# Per population comparison
ggplot(pi2, aes(x=Var1, y=Freq, fill=Var1))+
    geom_bar(stat="identity",position=position_dodge(width=0.8))+
    scale_fill_manual(values=div1)+
    theme_classic()+
    theme(axis.text.x = element_text(angle=45, hjust=1), legend.title = element_blank())+
    ggtitle(expression(paste("Changes in ", pi)))+
    xlab('')+ylab('Number of regions with P<0.05')
ggsave("../Output/Pi/Shuffle/Pi_significant_perComparison.png", width = 5, height = 4, dpi=300) 

#ordered
pi3<-pi2[order(pi2$Freq, decreasing = T),]
pi3$Var1<-factor(pi3$Var1, levels=paste0(unique(pi3$Var1)))
ggplot(pi3, aes(x=Var1, y=Freq, fill=Var1))+
    geom_bar(stat="identity",position=position_dodge(width=0.8))+
    scale_fill_manual(values=div2)+
    theme_classic()+
    theme(axis.text.x = element_text(angle=45, hjust=1), legend.title = element_blank())+
    ggtitle(expression(paste("Changes in ", pi, " (shuffle)")))+
    xlab('')+ylab('Number of regions with P<0.05')
ggsave("../Output/Pi/Shuffle/Pi_significant_perComparison_ordered.png", width = 5, height = 4, dpi=300) 

#Assign colors to the comparison group
names(div2)<- unique(pi3$Var1)

1.3.2 Theta

#Thetea
ths<-the[the$Var1 %in% yrs,]
ths$Var1<-factor(ths$Var1, levels=yrs)
ggplot(ths, aes(x=Var2, y=Freq, fill=Var1))+
    geom_bar(stat="identity",position=position_dodge(width=0.8))+
    scale_fill_manual(values=col3)+
    theme_classic()+
    theme(axis.text.x = element_text(angle=45, hjust=1), legend.title = element_blank())+
    ggtitle(paste0("Changes in theta"))+
    xlab('')+ylab('Number of regions with P<0.05')
ggsave("../Output/Pi/Shuffle/Theta_significant_perChrom_perPop.png", width = 5, height = 3, dpi=300) 

ggplot(the2, aes(x=Var1, y=Freq, fill=Var1))+
    geom_bar(stat="identity",position=position_dodge(width=0.8))+
    scale_fill_manual(values=div1)+
    theme_classic()+
    theme(axis.text.x = element_text(angle=45, hjust=1), legend.title = element_blank())+
    ggtitle(paste("Changes in theta"))+
    xlab('')+ylab('Number of regions with P<0.05')
ggsave("../Output/Pi/Shuffle/Theta_significant_perComparison.png", width = 5, height = 4, dpi=300) 


the3<-the2[order(the2$Freq, decreasing = T),]
the3$Var1<-factor(the3$Var1, levels=paste0(unique(the3$Var1)))

ggplot(the3, aes(x=Var1, y=Freq, fill=Var1))+
    geom_bar(stat="identity",position=position_dodge(width=0.8))+
    scale_fill_manual(name = "Var1",values = div2)+
    theme_classic()+
    theme(axis.text.x = element_text(angle=45, hjust=1), legend.title = element_blank())+
    ggtitle(paste("Changes in theta"))+
    xlab('')+ylab('Number of regions with P<0.05')
ggsave("../Output/Pi/Shuffle/Theta_significant_perComparison_ordered.png", width = 5, height = 4, dpi=300) 

1.3.3 Tajima’s D - No regions with P < 0.05

# Tajima's D
ggplot(Ds2, aes(x=Var1, y=Freq, fill=Var1))+
    geom_bar(stat="identity",position=position_dodge(width=0.8))+
    scale_fill_manual(values=div1)+
    theme_classic()+
    theme(axis.text.x = element_text(angle=45, hjust=1), legend.title = element_blank())+
    ggtitle(paste0("Changes in Tajima's D"))+
    xlab('')+ylab('Number of regions with P<0.05')
ggsave("../Output/Pi/Shuffle/TajimaD_significant_perComparison.png", width = 5, height = 4, dpi=300) 


D2<-D[D$Var1 %in% yrs,]
D2$Var1<-factor(D2$Var1, levels=yrs)
ggplot(D2, aes(x=Var2, y=Freq, fill=Var1))+
    geom_bar(stat="identity",position=position_dodge(width=0.8))+
    scale_fill_manual(values=col3)+
    theme_classic()+
    theme(axis.text.x = element_text(angle=45, hjust=1), legend.title = element_blank())+
     ggtitle(paste0("Changes in Tajima's D"))+
    xlab('')+ylab('Number of regions with P>0.05')
#ggsave("../Output/Pi/Shuffle/TajimaD_significant_perChrom_perPop.png", width = 8, height = 4, dpi=300) 
  • Most differences exist between 1996 and 2007
  • Chr25 has the most significant regions for changes in Pi and Theta



2 Fst genome-scan using a null model

  • from Pinsky et al. 2021 https://github.com/pinskylab/codEvol
  • Create a null distribution of expected maximum Fst
  • Define a genome-wide p-value for each region as (r+1)/(n+1), where r was the number of null simulations with maximum change greater than or equal to the empirical value and n was the number of shuffles

2.1 Shuffle Fst and create a null model of max Fst

  • Based on angsd persite Fst output
  1. After running angsd (realSFS), print a persite Fst file
```bash
#run FstPWSprint.sh

/home/jamcgirr/apps/angsd/misc/realSFS fst print  /home/ktist/ph/data/angsd/analysis/fst_PWS91_PWS07_persite_maf00.fst.idx > /home/ktist/ph/data/angsd/analysis/fst_PWS91_PWS07_persite_maf00.txt
bgzip /home/ktist/ph/data/angsd/analysis/fst_PWS91_PWS07_persite_maf00.txt

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


2. Run Fst_shuffle_pws.R at Farm (slurm script: angsd_fst_siteshuffle_null.sh)  


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjpbIiMgRnJvbSBodHRwczovL2dpdGh1Yi5jb20vcGluc2t5bGFiL2NvZEV2b2wgYW5nc2RfZnN0X3NpdGVzaHVmZmxlX251bGwuciIsIiIsIiMgc2h1ZmZsZSBBTkdTRCBwZXJzaXRlIEZTVCBBIGFuZCBCIHZhbHVlcyBhY3Jvc3Mgc2l0ZXMgYW5kIGNhbGN1bGF0ZSB3aW5kb3dlZCBGU1QgdG8gZ2V0IGEgbnVsbCBkaXN0cmlidXRpb24gb2YgbWF4IGdlbm9tZS13aWRlIEZTVCIsIiMgbmVlZCBmc3QgcGVyc2l0ZSBvdXRwdXQgZnJvbSBhbmdzZCBmc3RfKl9wZXJzaXRlX21hZjAwLnR4dC5neiBmaWxlcyIsIiIsIiMgdG8gcnVuIG9uIGZhcm0iLCIiLCIiLCIjIyMjIyBSIGNvZGUgICMjIyMjIyMgIiwiIiwiIyBwYXJhbWV0ZXJzIiwid2luc3ogPC0gNTAwMDAgIyB3aW5kb3cgc2l6ZSIsIndpbnN0cCA8LSAxMDAwMCAjIHdpbmRvdyBzdGVwIiwibnJlcCA8LSAxMDAwICMgbnVtYmVyIG9mIHJlc2h1ZmZsZXMiLCJtaW5sb2NpIDwtIDEwICMgbWluaW11bSBudW1iZXIgb2YgbG9jaSBwZXIgd2luZG93IHRvIGNvbnNpZGVyIiwiIiwiIiwiIyBsb2FkIGZ1bmN0aW9ucyIsInJlcXVpcmUoZGF0YS50YWJsZSkiLCIiLCIiLCIjIyMjIyMjIyMjIyMjIiwiIyBQcmVwIGRhdGEiLCIjIyMjIyMjIyMjIyMjIiwiIyBsb2FkIGZzdCBBL0IgZGF0YSIsIiNjYW4gPC0gZnJlYWQoJ2FuYWx5c2lzL0Nhbl80MC5DYW5fMTQuZnN0LkFCLmd6JykiLCIjc2V0bmFtZXMoY2FuLCBjKCdDSFJPTScsICdQT1MnLCAnQScsICdCJykpIiwiIiwicHdzOTE5NiA8LSBmcmVhZCgnL2hvbWUva3Rpc3QvcGgvZGF0YS9hbmdzZC9hbmFseXNpcy9mc3RfUFdTOTFfUFdTOTZfcGVyc2l0ZV9tYWYwMC50eHQuZ3onKSIsInNldG5hbWVzKHB3czkxOTYsIGMoJ0NIUk9NJywgJ1BPUycsICdBJywgJ0InKSkiLCJwd3M5MTA3IDwtIGZyZWFkKCcvaG9tZS9rdGlzdC9waC9kYXRhL2FuZ3NkL2FuYWx5c2lzL2ZzdF9QV1M5MV9QV1MwN19wZXJzaXRlX21hZjAwLnR4dC5neicpIiwic2V0bmFtZXMocHdzOTEwNywgYygnQ0hST00nLCAnUE9TJywgJ0EnLCAnQicpKSIsInB3czkxMTcgPC0gZnJlYWQoJy9ob21lL2t0aXN0L3BoL2RhdGEvYW5nc2QvYW5hbHlzaXMvZnN0X1BXUzkxX1BXUzE3X3BlcnNpdGVfbWFmMDAudHh0Lmd6JykiLCJzZXRuYW1lcyhwd3M5MTE3LCBjKCdDSFJPTScsICdQT1MnLCAnQScsICdCJykpIiwicHdzOTYwNyA8LSBmcmVhZCgnL2hvbWUva3Rpc3QvcGgvZGF0YS9hbmdzZC9hbmFseXNpcy9mc3RfUFdTOTZfUFdTMDdfcGVyc2l0ZV9tYWYwMC50eHQuZ3onKSIsInNldG5hbWVzKHB3czk2MDcsIGMoJ0NIUk9NJywgJ1BPUycsICdBJywgJ0InKSkiLCJwd3M5NjE3IDwtIGZyZWFkKCcvaG9tZS9rdGlzdC9waC9kYXRhL2FuZ3NkL2FuYWx5c2lzL2ZzdF9QV1M5Nl9QV1MxN19wZXJzaXRlX21hZjAwLnR4dC5neicpIiwic2V0bmFtZXMocHdzOTYxNywgYygnQ0hST00nLCAnUE9TJywgJ0EnLCAnQicpKSIsInB3czA3MTcgPC0gZnJlYWQoJy9ob21lL2t0aXN0L3BoL2RhdGEvYW5nc2QvYW5hbHlzaXMvZnN0X1BXUzA3X1BXUzE3X3BlcnNpdGVfbWFmMDAudHh0Lmd6JykiLCJzZXRuYW1lcyhwd3MwNzE3LCBjKCdDSFJPTScsICdQT1MnLCAnQScsICdCJykpIiwiIiwiUFdTPC1saXN0KCkiLCJQV1NbWzFdXTwtcHdzOTE5NiIsIlBXU1tbMl1dPC1wd3M5MTA3IiwiUFdTW1szXV08LXB3czkxMTciLCJQV1NbWzRdXTwtcHdzOTYwNyIsIlBXU1tbNV1dPC1wd3M5NjE3ICIsIlBXU1tbNl1dPC1wd3MwNzE3ICIsIiAgICAiLCIgIiwiIyBjcmVhdGUgbmV3IGNvbHVtbnMgYXMgaW5kaWNlcyBmb3Igd2luZG93cyIsIiIsImZvciAoaSBpbiAxOiBsZW5ndGgoUFdTKSl7IiwiICAgIGRmPC1QV1NbW2ldXSIsIiAgICBmb3IoaiBpbiAxOih3aW5zei93aW5zdHApKXsiLCIgICAgICAgIGRmWywgKHBhc3RlMCgnd2luJywgaikpIDo9IGZsb29yKChQT1MgLSAoai0xKSp3aW5zdHApL3dpbnN6KSp3aW5zeiArIHdpbnN6LzIgKyAoai0xKSp3aW5zdHBdIiwiICAgIH0iLCIgICAgUFdTW1tpXV08LWRmIiwifSIsIiIsIiMgbWFyayB3aW5kb3dzIHdpdGggPCBtaW5sb2NpIGZvciByZW1vdmFsIiwicmVtIDwtIHJlcCgwLCA2KSAjIG51bWJlciBvZiB3aW5kb3dzIHJlbW92ZWQgZm9yIGVhY2ggb2YgdGhlIDYgY29tcGFyaXNvbnMiLCIiLCJmb3IoaSBpbiAxOiBsZW5ndGgoUFdTKSl7IiwiICAgIHB3PC1QV1NbW2ldXSIsIiAgICBmb3IoaiBpbiAxOih3aW5zei93aW5zdHApKXsiLCIgICAgICAgIHB3d2luIDwtIHB3WywgLihuc25wcyA9IGxlbmd0aChQT1MpKSwgYnkgPSAuKHdpbiA9IGdldChwYXN0ZTAoJ3dpbicsIGopKSldICMgY2FsYyBudW0gc25wcyBwZXIgd2luZG93IiwiICAgICAgICByZW1bMV0gPC0gcmVtWzFdICsgcHd3aW5bLCBzdW0obnNucHMgPCBtaW5sb2NpKV0gIyByZWNvcmQgbnVtYmVyIHRvIGJlIHJlbW92ZWQiLCIgICAgICAgIHB3d2luWywgKHBhc3RlMCgnd2luJywgaiwgJ2tlZXAnKSkgOj0gMV0gIyBjcmVhdGUgY29sIHRvIG1hcmsgd2hpY2ggd2luZG93cyB0byBrZWVwIiwiICAgICAgICBwd3dpbltuc25wcyA8IG1pbmxvY2ksIChwYXN0ZTAoJ3dpbicsIGosICdrZWVwJykpIDo9IDBdICMgbWFyayB3aW5kb3dzIHRvIHJlbW92ZSIsIiAgICAgICAgcHd3aW5bLCBuc25wcyA6PSBOVUxMXSAjIGRyb3AgY29sdW1uIiwiICAgICAgICBzZXRuYW1lcyhwd3dpbiwgXCJ3aW5cIiwgcGFzdGUwKCd3aW4nLCBqKSkgIyBjaGFuZ2UgY29sIG5hbWUiLCIgICAgICAgIHB3IDwtIG1lcmdlKHB3LCBwd3dpbiwgYnkgPSBwYXN0ZTAoJ3dpbicsIGopLCBhbGwueCA9IFRSVUUpICMgbWVyZ2Uga2VlcGVyIGNvbCBiYWNrIHRvIGZ1bGwgZGF0YXNldCIsIiAgICB9IiwiICAgIFBXU1tbaV1dPC1wdyIsIn0iLCIiLCJyZW0gIyBudW1iZXIgb2Ygd2luZG93cyByZW1vdmVkIGZvciBlYWNoIGNvbXBhcmlzb24iLCIiLCIiLCIiLCIiLCIjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMiLCIjIHNodWZmbGUgYW5kIHJlY2FsYyB3aW5kb3dlZCBGU1QiLCIjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMiLCJjb2xubXMgPC0gYygnQ0hST00nLCAnUE9TJywgcGFzdGUwKCd3aW4nLCAxOih3aW5zei93aW5zdHApKSwgcGFzdGUwKCd3aW4nLCAxOih3aW5zei93aW5zdHApLCAna2VlcCcpKSAjIGxpc3Qgb2YgY29sdW1uIG5hbWVzIHdlIHdhbnQgb3V0IG9mIHRoZSBiYXNlIGRhdGEudGFibGUiLCIiLCIjIFBXUzA3LTE3IiwicG9wczwtYyhcIlBXUzkxLjk2XCIsXCJQV1M5MS4wN1wiLFwiUFdTOTEuMTdcIixcIlBXUzk2LjA3XCIsXCJQV1M5Ni4xN1wiLFwiUFdTMDcuMTdcIikiLCIiLCIiLCJmb3IocCBpbiAxOmxlbmd0aChQV1MpKXsiLCIgICAgcHJpbnQocGFzdGUwKCdTdGFydGluZyAnLCBwb3BzW3BdKSkiLCIiLCIgICAgcHc8LVBXU1tbcF1dIiwiICAgIGZvcihpIGluIDE6bnJlcCl7IiwiICAgICAgICBjYXQoaSk7IGNhdCgnICcpIiwiICAgICAgICAjIGNyZWF0ZSBuZXcgZGF0YXNldCIsIiAgICAgICAgaW5kcyA8LSBzYW1wbGUoMTpucm93KHB3KSwgbnJvdyhwdyksIHJlcGxhY2UgPSBGQUxTRSkiLCIgICAgICAgIHRlbXAgPC0gY2JpbmQocHdbLCAuLmNvbG5tc10sIHB3W2luZHMsIC4oQSwgQildKSAjIHNodWZmbGUgRlNUcyBhY3Jvc3MgcG9zaXRpb25zIiwiICAgICAgICAiLCIgICAgICAgICMgY2FsYyBmc3QgZm9yIGVhY2ggd2luZG93IHRvIGtlZXAiLCIgICAgICAgIGZvcihqIGluIDE6KHdpbnN6L3dpbnN0cCkpeyIsIiAgICAgICAgICAgIHRlbXAyIDwtIHRlbXBbZ2V0KHBhc3RlMCgnd2luJywgaiwgJ2tlZXAnKSkgPT0gMSwgXSAjIHRyaW0gdG8gd2luZG93cyB0byBrZWVwLiBjYW4ndCBjb21iaW5lIHdpdGggbmV4dCBsaW5lIGZvciBzb21lIHJlYXNvbi4iLCIgICAgICAgICAgICBpZihqID09MSkgdGVtcGZzdHMgPC0gdGVtcDJbLCAuKGZzdCA9IHN1bShBKS9zdW0oQikpLCBieSA9IC4oQ0hST00sIFBPUyA9IGdldChwYXN0ZTAoJ3dpbicsIGopKSldIiwiICAgICAgICAgICAgaWYoaiA+IDEpIHRlbXBmc3RzIDwtIHJiaW5kKHRlbXBmc3RzLCB0ZW1wMlssIC4oZnN0ID0gc3VtKEEpL3N1bShCKSksIGJ5ID0gLihDSFJPTSwgUE9TID0gZ2V0KHBhc3RlMCgnd2luJywgaikpKV0pIiwiICAgICAgICB9IiwiICAgICAgICAiLCIgICAgICAgICMgc2F2ZSB0aGUgbWF4IHdpbmRvd2VkIGZzdCIsIiAgICAgICAgIyBleGNsdWRlIHdpbmRvd3Mgd2l0aCBuZWdhdGl2ZSBtaWRwb2ludHMiLCIgICAgICAgIGlmKGkgPT0gMSkgbWF4ZnN0IDwtIHRlbXBmc3RzW1BPUyA+IDAsIG1heChmc3QsIG5hLnJtID0gVFJVRSldXHQiLCIgICAgICAgIGlmKGkgPiAxKSBtYXhmc3QgPC0gYyhtYXhmc3QsIHRlbXBmc3RzW1BPUyA+IDAsIG1heChmc3QsIG5hLnJtID0gVFJVRSldKSIsIiAgICB9IiwiIiwiICAgIHByaW50KHBhc3RlKCdNYXg6JywgbWF4KG1heGZzdCwgbmEucm0gPSBUUlVFKSwgJzsgOTV0aDonLCBxdWFudGlsZShtYXhmc3QsIHByb2IgPSAwLjk1LCBuYS5ybSA9IFRSVUUpKSkiLCIgICAgIiwiICAgIHdyaXRlLmNzdihtYXhmc3QsIGd6ZmlsZShwYXN0ZTAoJy9ob21lL2t0aXN0L3BoL2RhdGEvYW5nc2QvYW5hbHlzaXMvJyxwb3BzW3BdLCAnX3NpdGVzaHV1ZmZsZS5jc3YuZ3onKSksIHJvdy5uYW1lcyA9IEZBTFNFKSIsIiAgICBybShtYXhmc3QpIiwifSJdfQ== -->

```r
# From https://github.com/pinskylab/codEvol angsd_fst_siteshuffle_null.r

# shuffle ANGSD persite FST A and B values across sites and calculate windowed FST to get a null distribution of max genome-wide FST
# need fst persite output from angsd fst_*_persite_maf00.txt.gz files

# to run on farm


##### R code  ####### 

# parameters
winsz <- 50000 # window size
winstp <- 10000 # window step
nrep <- 1000 # number of reshuffles
minloci <- 10 # minimum number of loci per window to consider


# load functions
require(data.table)


#############
# Prep data
#############
# load fst A/B data
#can <- fread('analysis/Can_40.Can_14.fst.AB.gz')
#setnames(can, c('CHROM', 'POS', 'A', 'B'))

pws9196 <- fread('/home/ktist/ph/data/angsd/analysis/fst_PWS91_PWS96_persite_maf00.txt.gz')
setnames(pws9196, c('CHROM', 'POS', 'A', 'B'))
pws9107 <- fread('/home/ktist/ph/data/angsd/analysis/fst_PWS91_PWS07_persite_maf00.txt.gz')
setnames(pws9107, c('CHROM', 'POS', 'A', 'B'))
pws9117 <- fread('/home/ktist/ph/data/angsd/analysis/fst_PWS91_PWS17_persite_maf00.txt.gz')
setnames(pws9117, c('CHROM', 'POS', 'A', 'B'))
pws9607 <- fread('/home/ktist/ph/data/angsd/analysis/fst_PWS96_PWS07_persite_maf00.txt.gz')
setnames(pws9607, c('CHROM', 'POS', 'A', 'B'))
pws9617 <- fread('/home/ktist/ph/data/angsd/analysis/fst_PWS96_PWS17_persite_maf00.txt.gz')
setnames(pws9617, c('CHROM', 'POS', 'A', 'B'))
pws0717 <- fread('/home/ktist/ph/data/angsd/analysis/fst_PWS07_PWS17_persite_maf00.txt.gz')
setnames(pws0717, c('CHROM', 'POS', 'A', 'B'))

PWS<-list()
PWS[[1]]<-pws9196
PWS[[2]]<-pws9107
PWS[[3]]<-pws9117
PWS[[4]]<-pws9607
PWS[[5]]<-pws9617 
PWS[[6]]<-pws0717 
    
 
# create new columns as indices for windows

for (i in 1: length(PWS)){
    df<-PWS[[i]]
    for(j in 1:(winsz/winstp)){
        df[, (paste0('win', j)) := floor((POS - (j-1)*winstp)/winsz)*winsz + winsz/2 + (j-1)*winstp]
    }
    PWS[[i]]<-df
}

# mark windows with < minloci for removal
rem <- rep(0, 6) # number of windows removed for each of the 6 comparisons

for(i in 1: length(PWS)){
    pw<-PWS[[i]]
    for(j in 1:(winsz/winstp)){
        pwwin <- pw[, .(nsnps = length(POS)), by = .(win = get(paste0('win', j)))] # calc num snps per window
        rem[1] <- rem[1] + pwwin[, sum(nsnps < minloci)] # record number to be removed
        pwwin[, (paste0('win', j, 'keep')) := 1] # create col to mark which windows to keep
        pwwin[nsnps < minloci, (paste0('win', j, 'keep')) := 0] # mark windows to remove
        pwwin[, nsnps := NULL] # drop column
        setnames(pwwin, "win", paste0('win', j)) # change col name
        pw <- merge(pw, pwwin, by = paste0('win', j), all.x = TRUE) # merge keeper col back to full dataset
    }
    PWS[[i]]<-pw
}

rem # number of windows removed for each comparison




####################################
# shuffle and recalc windowed FST
####################################
colnms <- c('CHROM', 'POS', paste0('win', 1:(winsz/winstp)), paste0('win', 1:(winsz/winstp), 'keep')) # list of column names we want out of the base data.table

# PWS07-17
pops<-c("PWS91.96","PWS91.07","PWS91.17","PWS96.07","PWS96.17","PWS07.17")


for(p in 1:length(PWS)){
    print(paste0('Starting ', pops[p]))

    pw<-PWS[[p]]
    for(i in 1:nrep){
        cat(i); cat(' ')
        # create new dataset
        inds <- sample(1:nrow(pw), nrow(pw), replace = FALSE)
        temp <- cbind(pw[, ..colnms], pw[inds, .(A, B)]) # shuffle FSTs across positions
        
        # calc fst for each window to keep
        for(j in 1:(winsz/winstp)){
            temp2 <- temp[get(paste0('win', j, 'keep')) == 1, ] # trim to windows to keep. can't combine with next line for some reason.
            if(j ==1) tempfsts <- temp2[, .(fst = sum(A)/sum(B)), by = .(CHROM, POS = get(paste0('win', j)))]
            if(j > 1) tempfsts <- rbind(tempfsts, temp2[, .(fst = sum(A)/sum(B)), by = .(CHROM, POS = get(paste0('win', j)))])
        }
        
        # save the max windowed fst
        # exclude windows with negative midpoints
        if(i == 1) maxfst <- tempfsts[POS > 0, max(fst, na.rm = TRUE)]  
        if(i > 1) maxfst <- c(maxfst, tempfsts[POS > 0, max(fst, na.rm = TRUE)])
    }

    print(paste('Max:', max(maxfst, na.rm = TRUE), '; 95th:', quantile(maxfst, prob = 0.95, na.rm = TRUE)))
    
    write.csv(maxfst, gzfile(paste0('/home/ktist/ph/data/angsd/analysis/',pops[p], '_siteshuuffle.csv.gz')), row.names = FALSE)
    rm(maxfst)
}

2.2 Assess the results from site reshuffled Fst

Need - sitehuffle results fst_siteshuuffle.csv.gz - window-based Fst results from angsd _50kWindow_maf00 - persite Fst files after LD pruning (with AB info) *_persite_maf00_ldtrim.csv.gz

#from https://github.com/pinskylab/codEvol

# parameters
minloci <- 10 
winsz <- 50000 # window size
winstp <- 10000 # window step

# load functions
require(data.table)
require(ggplot2)

calcp <- function(fst, null) return((sum(null > fst)+1)/(length(null)+1)) # equation from North et al. 2002 Am J Hum Gen

# read in and prep data
years<-c("PWS91","PWS96","PWS07","PWS17")
comb<-t(combn(years, 2))

# continuous nucleotide position for the whole genome
chrmax <- fread('../Data/new_vcf/chr_sizes.bed')
chrmax<-chrmax[,-2]
colnames(chrmax)<-c("chr", "len")
chrmax$start<-c(0,cumsum(chrmax$len)[1:(nrow(chrmax)-1)])
setkey(chrmax, chr)

prunedFst_gw<-data.frame(pop1=comb[,1],pop2=comb[,2])
#for (p in 1: nrow(comb)){
    for (p in 1:5){
    #genome-wide max Fst from reshuffling (unlinked sites)
    null<-fread(paste0('Data/shuffle/', comb[p,1],".",comb[p,2],'_fst_siteshuuffle.csv.gz'))
    #window-based Fst from angsd            
    fstwin <- fread(paste0('Data/new_vcf/angsd/fromVCF/2D/fst_',comb[p,1],"_",comb[p,2],'_50kWindow_maf00'), header = FALSE, col.names = c('region', 'chr', 'midPos', 'Nsites', 'fst')) # all sites
    #persite fst (AB) -pruned sites
    AB <- fread(paste0('Data/shuffle/fst_',comb[p,1],"_",comb[p,2],'_persite_maf00_ldtrim.csv.gz'))
    
    # merge nucleotide position into the frequency files
    setkey(AB, CHROM)
    AB <- AB[chrmax[, .(CHROM = chr, start)], ]
    AB[, posgen := POS + start]
    AB[,start := NULL]
    
    # Calc genome-wide FST
    prunedFst_gw$gwFst[p]<-AB[!is.na(A), sum(A)/sum(B)]
    
    # Calc windowed FST
    # create new columns as indices for windows
    for(j in 1:(winsz/winstp)){
        AB[, (paste0('win', j)) := floor((POS - (j-1)*winstp)/winsz)*winsz + winsz/2 + (j-1)*winstp]
    }
    
    # calc fst and # snps per window
    for(j in 1:(winsz/winstp)){
        if(j ==1){
            fstwin <- AB[, .(fst = sum(A)/sum(B), nloci = length(POS)), by = .(CHROM, midPos = get(paste0('win', j)))]
        } 
        if(j > 1){
            fstwin <- rbind(fstwin, AB[, .(fst = sum(A)/sum(B), nloci = length(POS)), by = .(CHROM, midPos = get(paste0('win', j)))])
        } 
    }
    
    ## null model stats
    fstats<-null[, .(max = max(x), u95 = quantile(x, probs = 0.95))]
    prunedFst_gw$null.Fstmax[p]<-fstats$max[1]
    prunedFst_gw$null.Fst.95q[p]<-fstats$u95[1]
    
    #Calculate p-values
    pval<-fstwin[, p := calcp(fst, null$x), by = .(CHROM, midPos)]
    
    #combined the datasets
    pair<-paste0(comb[p,1],"_",comb[p,2])
    fstwin[, pop := pair ]
    write.csv(fstwin, paste0("../Output/Fst/fst_siteshuffle.pairwise_", pair,".csv"))
}

write.csv(prunedFst_gw, "../Output/Fst/fst_siteshuffle_pws_summary.csv")

data<-data.frame()
for(i in 1: nrow(comb)){
    pair<-paste0(comb[i,1],"_",comb[i,2])
    df<-fread(paste0("../Output/Fst/fst_siteshuffle.pairwise_", pair,".csv"))
    df<-df[!is.na(fst) & midPos > 0, ]
    df[,ch:= as.integer(gsub("chr",'', CHROM))]
    setkey(df,CHROM)
    df<-df[chrmax[, .(CHROM = chr, start)], ]
    df[,pos.gw:= midPos+start]
    #setkey(df, ch, midPos)
    data<-rbind(data, df)
}

write.csv(data, file = gzfile('../Output/Fst/fst_siteshuffle.angsd.PWS.csv.gz'), row.names = FALSE)

2.3 Plot the results

data<-fread('../Output/Fst/fst_siteshuffle.angsd.PWS.csv.gz')

#pop as factor 
data[, pop := factor(pop, levels = c("PWS91_PWS96","PWS91_PWS07","PWS91_PWS17","PWS96_PWS07","PWS96_PWS17","PWS07_PWS17"))]

# plot p-value vs. position
two_cols<-rep(c("steelblue","lightblue"), times=13)
minloci <- 10 # minimum number of loci per window to consider

ggplot(data[nloci >= minloci, ], aes(pos.gw, -log10(p), color = factor(ch))) + 
    geom_point(size = 0.2, alpha = 0.5) +
    facet_wrap(~pop, ncol = 1) +
    scale_color_manual(values = two_cols, guide='none') +
    geom_hline(yintercept = -log10(0.05), linetype = 'dashed', color = 'grey')+
    theme_bw()+
    theme(axis.text.x = element_blank())+
    xlab("Genome")
ggsave("../Output/Fst/fst.siteshuffle.p_vs_pos.PWS.png", width = 7.5, height = 7.5,dpi = 300)

  • Number of loci per window vs. Fst p-values
# plot p-value vs. nloci to check if any biases exist
ggplot(data, aes(nloci, -log10(p), color = pop)) + 
    geom_point(size = 0.6, alpha = 0.5) +
    geom_hline(yintercept = -log10(0.05), linetype = 'dashed', color = 'grey') + 
    scale_x_log10()+theme_bw()+theme(legend.title = element_blank())+
    guides(color = guide_legend(override.aes = list(size = 3, alpha=1)))
ggsave("../Output/Fst/fst.siteshuffle.p_vs_nloci.PWS.png", width = 8, height =4,dpi = 300)

2.4 Significant outlier regions

2.4.1 p-value < 0.05

```r
outliers<-data[p < 0.05, ]
write.csv(outliers, \../Output/Fst/Fst_nullshuffling_outliers_PWS.csv\)

counts<-data.table(table(outliers$pop))
ggplot(counts, aes(x=V1, y=N))+
    geom_bar(stat=\identity\, fill=cols[4])+
    xlab('')+ylab('')+
    theme_classic()+
    theme(axis.text.x = element_text(angle=45, hjust=1))
ggsave(\../Output/Fst/Fst_outliers_byComparison.png\, width = 4, height = 3, dpi=300 )


counts$V1<-factor(counts$V1, levels=c(\PWS96_PWS07\,\PWS07_PWS17\,\PWS91_PWS96\, \PWS91_PWS07\, \PWS91_PWS17\,\PWS96_PWS17\))
ggplot(counts, aes(x=V1, y=N, fill=V1))+
    geom_bar(stat=\identity\)+
    scale_fill_manual(name = \V\,values = div2)+
    xlab('')+ylab('Number of regions with P<0.05')+
    theme_classic()+ggtitle(\Fst (shuffle)\)+
    theme(axis.text.x = element_text(angle=45, hjust=1), legend.title = element_blank())
ggsave(\../Output/Fst/Fst_outliers_byComparison_ordered.png\, width = 4, height = 3, dpi=300 )



#outliers<-read.csv(\../Output/Fst/Fst_nullshuffling_outliers_PWS.csv\, row.names = 1)
knitr::kable(outliers[,c(2,3,4,5,6,7,8)], caption=\Outlier regions\)

#Outlier regions
#CHROM  midPos  fst nloci   p   pop ch
#chr25  825000  0.1210997   138 0.003996    PWS96_PWS07 25
#chr25  835000  0.0818802   279 0.020979    PWS96_PWS07 25
#chr25  845000  0.0721926   311 0.033966    PWS96_PWS07 25
#chr25  805000  0.1242652   166 0.002997    PWS96_PWS07 25
#chr25  815000  0.1377218   155 0.002997    PWS96_PWS07 25
#chr19  1375000 0.1251602   42  0.003996    PWS07_PWS17 19

```

2.4.2 Top 1% outliers & P< 0.2

  • relaxing the cutoff


* NO overlapping regions from PCAngsd selection scan

LS0tCnRpdGxlOiAiU2h1ZmZsaW5nIFBpL0ZzdC9UaGV0YS9EIgpvdXRwdXQ6CiAgaHRtbF9ub3RlYm9vazoKICAgICAgdG9jOiB0cnVlIAogICAgICB0b2NfZmxvYXQ6IHRydWUKICAgICAgbnVtYmVyX3NlY3Rpb25zOiBUcnVlCiAgICAgIHRoZW1lOiBsdW1lbgogICAgICBoaWdobGlnaHQ6IHRhbmdvCiAgICAgIGNvZGVfZm9sZGluZzogaGlkZQogICAgICBkZl9wcmludDogcGFnZWQKLS0tCgojIEZpbmQgc2lnbmlmaWNhbnRseSBkaWZmZXJlbnQgbG9jaS9yZWdpb25zIGJldHdlZW4geWVhcnMgdXNpbmcgc2h1ZmZsaW5nIHRvIGNyZWF0ZSBudWxsIGRpc3RyaWJ1dGlvbnN7LnVubnVtYmVyZWR9ICAKIApgYGB7ciBldmFsPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpzb3VyY2UoIi4uL1JzY3JpcHRzL0Jhc2VTY3JpcHRzLlIiKQpyZXF1aXJlKGRhdGEudGFibGUpCnJlcXVpcmUocGx5cikKcmVxdWlyZShSQ29sb3JCcmV3ZXIpCmBgYAoKCiMgVGhldGEgbnVsbCBkaXN0cmlidXRpb24KCiMjIENyZWF0ZSBwZXJzaXRlIHRoZXRhLnBlc3QuZ3ogZmlsZXMgaW4gYW5nc2QKCmBgYHtiYXNofQoKL2hvbWUvamFtY2dpcnIvYXBwcy9hbmdzZC9taXNjL3RoZXRhU3RhdCAgcHJpbnQgL2hvbWUva3Rpc3QvcGgvZGF0YS9hbmdzZC90aGV0YS9QV1M5MV9tYWYwMC50aGV0YXMuaWR4IHwgZ3ppcCA+IC9ob21lL2t0aXN0L3BoL2RhdGEvYW5nc2QvdGhldGEvUFdTOTFfbWFmMDBfdGhldGFzLnBlc3RQRy5negoKIyBSdW4gcHJpbnRUaGV0YS5zaCBhdCBmYXJtCmBgYAoKCiMjIFJ1biBSIHNjcmlwdHMgYXQgZmFybSB0byBjcmVhdGUgbnVsbCBkaXN0cmlidXRpb24gb2YgZ2Vub21lLXdpZGUgdGhldGFzCgpSdW4gYW5nc2RfdGhldGFfc2l0ZXNodWZmbGVfbnVsbC5zaCBhdCBmYXJtLCB3aGljaCBydW5zIFBpX3NodWZmbGVfcHdzLlIgCiAtIFRha2VzIGxvbmcgdGltZSwgc28gY3JlYXRlIGEgc2NyaXB0IGZvciBlYWNoIHBvcC55ZWFyIGNvbWJpbmF0aW9uIGFuZCBydW4gc2VwYXJhdGVseSAKCk91dHB1dCBmcm9tIHRoZXRhIHNodWZmbGluZyByZXN1bHRzIGFyZSBpbiBEYXRhL3NodWZmbGUvdGhldGEuc2l0ZXNodWZmbGUuNTAwMDAuUFdTOTFfUFdTOTYuY3N2Lmd6CgoKIyMgQ2FsY3VsYXRlIHAtdmFsdWVzIGZvciBkaWZmZXJlbmNlcyBpbiB0aGV0YSwgcGksIGFuZCBUYWppbWEncyBEIHVzaW5nIHRoZSByZXN1bHRzIGZyb20gc2h1ZmZsaW5nIAoKYGBge3IgZXZhbD1GQUxTRSwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KI2Zyb20gY29kRXZvbCBhbmdzZF90aGV0YV9zaXRlc2h1ZmZsZV9udWxsX3N0YXRzLlIKCiMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIwojIGxvYWQgZnVuY3Rpb25zCiMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIwpyZXF1aXJlKGRhdGEudGFibGUpCnJlcXVpcmUocGx5cikKcmVxdWlyZShnZ3Bsb3QyKQpyZXF1aXJlKFJDb2xvckJyZXdlcikKCgojIyMjIyMjIyMjIyMjIyMjIyMjIyMKIyByZWFkIGluIGFuZCBwcmVwIGRhdGEKIyMjIyMjIyMjIyMjIyMjIyMjIyMjCgp5ZWFyczwtYygiUFdTOTEiLCJQV1M5NiIsIlBXUzA3IiwiUFdTMTciKQpjb21iPC10KGNvbWJuKHllYXJzLCAyKSkKCmNocm1heCA8LSBmcmVhZCgnLi4vRGF0YS9uZXdfdmNmL2Nocl9zaXplcy5iZWQnKQpjaHJtYXg8LWNocm1heFssLTJdCmNvbG5hbWVzKGNocm1heCk8LWMoImNociIsICJsZW4iKQpjaHJtYXgkc3RhcnQ8LWMoMCxjdW1zdW0oY2hybWF4JGxlbilbMToobnJvdyhjaHJtYXgpLTEpXSkKCmNocm1heCRlbmQ8LWN1bXN1bShjaHJtYXgkbGVuKQpjaHJtYXgkbWlkZGxlPC0oY2hybWF4JGVuZC1jaHJtYXgkc3RhcnQpLzIrY2hybWF4JHN0YXJ0Cgojc2V0a2V5KGNocm1heCwgY2hyKQoKI0Z1bmN0aW9ucyB0byBjYWxjdWxhdGUgcC12YWx1ZXMgZnJvbSBjb2RFdm9sCmNhbGNwRyA8LSBmdW5jdGlvbih0aGV0YWNoYW5nZSwgbnVsbCl7ICMgZm9yIGluY3JlYXNlcyBpbiB0aGV0YQogIHJldHVybigoc3VtKG51bGwgPiB0aGV0YWNoYW5nZSkrMSkvKGxlbmd0aChudWxsKSsxKSkgIyBlcXVhdGlvbiBmcm9tIE5vcnRoIGV0IGFsLiAyMDAyIEFtIEogSHVtIEdlbgp9CmNhbGNwTCA8LSBmdW5jdGlvbih0aGV0YWNoYW5nZSwgbnVsbCl7ICMgZm9yIGRlY3JlYXNlcyBpbiB0aGV0YQogIHJldHVybigoc3VtKG51bGwgPCB0aGV0YWNoYW5nZSkrMSkvKGxlbmd0aChudWxsKSsxKSkgIyBlcXVhdGlvbiBmcm9tIE5vcnRoIGV0IGFsLiAyMDAyIEFtIEogSHVtIEdlbgp9CgoKY2hfY29scyA8LSBicmV3ZXIucGFsKDQsICdQYWlyZWQnKVtyZXAoMToyLDEzKV0KRGF0YWxsPC1kYXRhLnRhYmxlKCkKZm9yIChwIGluIDE6IG5yb3coY29tYikpewogICAgIyBtYXggdGhldGEgcGVyIGdlbm9tZSBmcm9tIHJlc2h1ZmZsaW5nIChhbGwgc2l0ZXMpIGZyb20gYW5nc2RfdGhldGFfc2l0ZXNodWZmbGVfbnVsbC5yCiAgICBudWxsPC1mcmVhZChwYXN0ZTAoJy4uL0RhdGEvc2h1ZmZsZS90aGV0YS5zaXRlc2h1ZmZsZS41MDAwMC4nLCBjb21iW3AsMV0sIl8iLGNvbWJbcCwyXSwnLmNzdi5neicpKQogICAgCiAgICAjdXBwZXIgYW5kIGxvd2VyIDk1JQogICAgbnVsbFssIC4odFdkX2w5NSA9IHF1YW50aWxlKG1pbnRXZCwgMC4wNSksIHRXZF91OTUgPSBxdWFudGlsZShtYXh0V2QsIHByb2JzID0gMC45NSksCiAgICAgICAgICAgICAgICB0UGRfbDk1ID0gcXVhbnRpbGUobWludFBkLCAwLjA1KSwgdFBkX3U5NSA9IHF1YW50aWxlKG1heHRQZCwgcHJvYnMgPSAwLjk1KSwKICAgICAgICAgICAgICAgIHREZF9sOTUgPSBxdWFudGlsZShtaW50RGQsIDAuMDUpLCB0RGRfdTk1ID0gcXVhbnRpbGUobWF4dERkLCBwcm9icyA9IDAuOTUpKV0KCiAgICAjYXNzaWduKHBhc3RlMCgibnVsbC4iLGNvbWJbcCwxXSwiXyIsY29tYltwLDJdKSwgbnVsbCkgICAgCiAgICAKICAgICMgc2xpZGluZyB3aW5kb3dzIHRoZXRhIGNoYW5nZSAoR0FUSyBzaXRlcykgZnJvbSBhbmdzZF90aGV0YV9zaXRlc2h1ZmZsZV9udWxsLnIKICAgIGRhdDwtZnJlYWQocGFzdGUwKCcuLi9EYXRhL3NodWZmbGUvdGhldGFfY2hhbmdlX3JlZ2lvbl81MDAwMC4nLCBjb21iW3AsMV0sIl8iLGNvbWJbcCwyXSwnLmNzdi5neicpLCBkcm9wID0gMSkKICAgIGRhdFsscG9wOj1wYXN0ZTAoY29tYltwLDFdLCJfIixjb21iW3AsMl0pXQogICAgCiAgICBkYXQ8LW1lcmdlKGRhdCwgY2hybWF4WyxjKCJjaHIiLCJzdGFydCIpXSwgYnkueD0iQ2hyb21vIiwgYnkueSA9ICJjaHIiKQogICAgZGF0WywgUE9TZ2VuIDo9IFdpbkNlbnRlciArIHN0YXJ0XQogICAgZGF0WyxzdGFydCA6PSBOVUxMXSAjcmVtb3ZlIHN0YXJ0CiAgICAKICAgICNjYWxjdWxhdGUgcC12YWx1ZXMKICAgICMxLiB0aGV0YVcgbG9jaQogICAgZGF0W3RXZCA+IDAsIHRXZC5wIDo9IGNhbGNwRyh0V2QsIG51bGwkbWF4dFdkKSwgYnkgPSAuKENocm9tbywgV2luQ2VudGVyKV0gIyB0aGV0YVcgIGxvY2kKICAgIGRhdFt0V2QgPD0gMCwgdFdkLnAgOj0gY2FsY3BMKHRXZCwgbnVsbCRtaW50V2QpLCBieSA9IC4oQ2hyb21vLCBXaW5DZW50ZXIpXQogICAgIzIuIHRoZXRhIHBpCiAgICBkYXRbdFBkID4gMCwgdFBkLnAgOj0gY2FsY3BHKHRQZCwgbnVsbCRtYXh0UGQpLCBieSA9IC4oQ2hyb21vLCBXaW5DZW50ZXIpXSAjIHRoZXRhIHBpCiAgICBkYXRbdFBkIDw9IDAsIHRQZC5wIDo9IGNhbGNwTCh0UGQsIG51bGwkbWludFBkKSwgYnkgPSAuKENocm9tbywgV2luQ2VudGVyKV0KICAgIAogICAgI1RhamltYSdzIEQKICAgIGRhdFt0RGQgPiAwLCB0RGQucCA6PSBjYWxjcEcodERkLCBudWxsJG1heHREZCksIGJ5ID0gLihDaHJvbW8sIFdpbkNlbnRlcildICMgdGFqaW1hJ3MgRAogICAgZGF0W3REZCA8PSAwLCB0RGQucCA6PSBjYWxjcEwodERkLCBudWxsJG1pbnREZCksIGJ5ID0gLihDaHJvbW8sIFdpbkNlbnRlcildCgogICAgd3JpdGUuY3N2KGRhdCwgZmlsZT1nemZpbGUocGFzdGUwKCcuLi9PdXRwdXQvUGkvU2h1ZmZsZS90aGV0YV9zaXRlc2h1ZmZsZV8nLCBjb21iW3AsMV0sIl8iLGNvbWJbcCwyXSwnLmNzdi5neicpKSkKICAgIAogICAgRGF0YWxsPC1yYmluZChEYXRhbGwsIGRhdCkKfQoKd3JpdGUuY3N2KERhdGFsbCwgZmlsZT1nemZpbGUocGFzdGUwKCcuLi9PdXRwdXQvUGkvU2h1ZmZsZS90aGV0YV9zaXRlc2h1ZmZsZV9QV1Nfc3VtbWFyeS5jc3YuZ3onKSkpCmBgYCAgIAoKIyMgUGxvdCB0aGUgcmVzdWx0cyAgCgojIyMgQ2hhbmdlcyBpbiBQaS9UaGV0YS9EIGJldHdlZW4geWVhcnMgCmBgYHtyIGV2YWw9RkFMU0UsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CgojIyBQbG90IHRoZSByZXN1bHRzICMjCkRhdGFsbDwtZnJlYWQoJy4uL091dHB1dC9QaS9TaHVmZmxlL3RoZXRhX3NpdGVzaHVmZmxlX1BXU19zdW1tYXJ5LmNzdi5neicpCgp3aW5zeiA9IDVlNCAKCiNDaGFuZ2VzIGluIFBpIGJldHdlZW4geWVhcnMKRGF0YWxsJENocm9tbzwtZmFjdG9yKERhdGFsbCRDaHJvbW8sIGxldmVscz1jKHBhc3RlMCgiY2hyIiwxOjI2KSkpCkRhdGFsbCRwb3A8LWZhY3RvcihEYXRhbGwkcG9wLCBsZXZlbHM9YygiUFdTOTFfUFdTOTYiLCJQV1M5MV9QV1MwNyIsIlBXUzkxX1BXUzE3IiwiUFdTOTZfUFdTMDciLCJQV1M5Nl9QV1MxNyIsIlBXUzA3X1BXUzE3IikpCmNoX2NvbHMgPC0gYnJld2VyLnBhbCg0LCAnUGFpcmVkJylbcmVwKDE6MiwxMyldCmdncGxvdChEYXRhbGwsIGFlcyhQT1NnZW4sIHRQZC93aW5zeiwgY29sb3IgPSBDaHJvbW8pKSArIAogICAgZ2VvbV9wb2ludChzaXplID0gMC41LCBhbHBoYSA9IDAuMykgKwogICAgZmFjZXRfd3JhcCh+cG9wLCBuY29sID0gMSkgKwogICAgc2NhbGVfY29sb3JfbWFudWFsKHZhbHVlcyA9IGNoX2NvbHMsIGd1aWRlPSJub25lIikgKwogICAgeWxhYignQ2hhbmdlIGluIHBpIHBlciBzaXRlJykreGxhYigiQ2hyb21vc29tZSIpKwogICAgZ2d0aXRsZSgiQ2hhbmdlcyBpbiBQaSIpKwogICAgc2NhbGVfeF9jb250aW51b3VzKGJyZWFrcz1jaHJtYXgkbWlkZGxlLCBsYWJlbHM9MToyNikrCiAgICB0aGVtZV9idygpCmdnc2F2ZShwYXN0ZTAoJy4uL091dHB1dC9QaS9TaHVmZmxlL0NoYW5nZXNfaW5fUGlfUFdTLnBuZycpLCB3aWR0aCA9IDcuNSwgaGVpZ2h0ID0gOSwgZHBpID0gMzAwKQoKIyBwbG90IHRoZXRhX1dhdHRlcnNvbiBjaGFuZ2UKZ2dwbG90KERhdGFsbCwgYWVzKFBPU2dlbiwgdFdkL3dpbnN6LCBjb2xvciA9IENocm9tbykpICsgCiAgZ2VvbV9wb2ludChzaXplID0gMC41LCBhbHBoYSA9IDAuMykgKwogIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSBjaF9jb2xzLCBndWlkZT0ibm9uZSIpICsKICAgIGZhY2V0X3dyYXAofnBvcCwgbmNvbCA9IDEpICsKICB5bGFiKCdDaGFuZ2UgaW4gV2F0dGVyc29ucyB0aGV0YSBwZXIgc2l0ZScpK3hsYWIoIkNocm9tb3NvbWUiKSsKICBnZ3RpdGxlKCJDaGFuZ2VzIGluIFBpIikrCiAgICBzY2FsZV94X2NvbnRpbnVvdXMoYnJlYWtzPWNocm1heCRtaWRkbGUsIGxhYmVscz0xOjI2KSsKICAgIHRoZW1lX2J3KCkKZ2dzYXZlKCcuLi9PdXRwdXQvUGkvU2h1ZmZsZS9DaGFuZ2VzX2luX3RoZXRhV19QV1MucG5nJywgd2lkdGggPSA3LjUsIGhlaWdodCA9IDksIGRwaSA9IDMwMCkKCiMgcGxvdCBUYWppbWEncyBEIGNoYW5nZQpnZ3Bsb3QoRGF0YWxsLCBhZXMoUE9TZ2VuLCB0RGQsIGNvbG9yID0gQ2hyb21vKSkgKyAKICAgIGdlb21fcG9pbnQoc2l6ZSA9IDAuNSwgYWxwaGEgPSAwLjMpICsKICAgIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSBjaF9jb2xzLGd1aWRlPSJub25lIikgKwogICAgZmFjZXRfd3JhcCh+cG9wLCBuY29sID0gMSkgKwogICAgeWxhYignQ2hhbmdlIGluIFRhamltYXMgRCBwZXIgd2luZG93JykreGxhYigiQ2hyb21vc29tZSIpKwogICAgZ2d0aXRsZSgiQ2hhbmdlcyBpbiBUYWppbWEncyBEIikrCiAgICBzY2FsZV94X2NvbnRpbnVvdXMoYnJlYWtzPWNocm1heCRtaWRkbGUsIGxhYmVscz0xOjI2KSsKICAgIHRoZW1lX2J3KCkKZ2dzYXZlKCcuLi9PdXRwdXQvUGkvU2h1ZmZsZS9DaGFuZ2VzX2luX1RhamltYXNEX1BXUy5wbmcnLCB3aWR0aCA9IDcuNSwgaGVpZ2h0ID0gOSwgZHBpID0gMzAwKQoKYGBgCiFbXSguLi9PdXRwdXQvUGkvU2h1ZmZsZS9DaGFuZ2VzX2luX1BpX1BXUy5wbmcpe3dpZHRoPTcwJX0KCiFbXSguLi9PdXRwdXQvUGkvU2h1ZmZsZS9DaGFuZ2VzX2luX3RoZXRhV19QV1MucG5nKXt3aWR0aD03MCV9CgohW10oLi4vT3V0cHV0L1BpL1NodWZmbGUvQ2hhbmdlc19pbl9UYWppbWFzRF9QV1MucG5nKXt3aWR0aD03MCV9CgojIyMgUGxvdCB0aGUgcC12YWx1ZXMgZm9yIGRpZmZlcmVuY2VzIGluIFAvVGhldGEgYWxvbmcgdGhlIGdlbm9tZSAKCmBgYHtyIGV2YWw9RkFMU0UsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CiMgcGxvdCBwaSBwLXZhbHVlIHZzLiBwb3NpdGlvbiAKY2hfY29scyA8LSBicmV3ZXIucGFsKDQsICdQYWlyZWQnKVtyZXAoMToyLDEzKV0KRGF0YWxsJENocm9tbzwtZmFjdG9yKERhdGFsbCRDaHJvbW8sIGxldmVscz1jKHBhc3RlMCgiY2hyIiwxOjI2KSkpCkRhdGFsbCRwb3A8LWZhY3RvcihEYXRhbGwkcG9wLCBsZXZlbHM9YygiUFdTOTFfUFdTOTYiLCJQV1M5MV9QV1MwNyIsIlBXUzkxX1BXUzE3IiwiUFdTOTZfUFdTMDciLCJQV1M5Nl9QV1MxNyIsIlBXUzA3X1BXUzE3IikpCgpnZ3Bsb3QoRGF0YWxsLCBhZXMoUE9TZ2VuLCAtbG9nMTAodFBkLnApKnNpZ24odFBkKSwgY29sb3IgPSBDaHJvbW8pKSArIAogICAgZ2VvbV9wb2ludChzaXplID0gMC40LCBhbHBoYSA9IDAuMykgKwogICAgZmFjZXRfd3JhcCh+cG9wLCBuY29sID0gMSkgKwogICAgc2NhbGVfY29sb3JfbWFudWFsKHZhbHVlcyA9IGNoX2NvbHMsIGd1aWRlPSJub25lIikgK3hsYWIoIkNocm9tb3NvbWUiKSsKICAgIHlsYWIoImxvZzEwKFAtdmFsdWUpIikrCiAgICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSBsb2cxMCgwLjA1KSwgbGluZXR5cGUgPSAnZGFzaGVkJywgY29sb3IgPSAncmVkJyxzaXplPTAuMykgKwogICAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gLWxvZzEwKDAuMDUpLCBsaW5ldHlwZSA9ICdkYXNoZWQnLCBjb2xvciA9ICdyZWQnLCBzaXplPTAuMykrCiAgICBnZ3RpdGxlKCJQLXZhbHVlcyBmb3IgY2hhbmdlcyBpbiBQaSIpKwogICAgc2NhbGVfeF9jb250aW51b3VzKGJyZWFrcz1jaHJtYXgkbWlkZGxlLCBsYWJlbHM9MToyNikrCiAgICB0aGVtZV9idygpCmdnc2F2ZSgnLi4vT3V0cHV0L1BpL1NodWZmbGUvQ2hhbmdlc19pbl9QaS5zaXRlc2h1ZmZsZS5wLXZhbHVlc19QV1MucG5nJywgd2lkdGggPSA3LjUsIGhlaWdodCA9MTEsIHVuaXRzID0gJ2luJywgZHBpID0gMzAwKQoKCiMgcGxvdCB0aGV0YVcgcC12YWx1ZSB2cy4gcG9zaXRpb24gKGFsbCBsb2NpKQpnZ3Bsb3QoRGF0YWxsLCBhZXMoUE9TZ2VuLCAtbG9nMTAodFdkLnApKnNpZ24odFdkKSwgY29sb3IgPSBDaHJvbW8pKSArIAogIGdlb21fcG9pbnQoc2l6ZSA9IDAuNCwgYWxwaGEgPSAwLjMpICsKICBmYWNldF93cmFwKH5wb3AsIG5jb2wgPSAxKSArCiAgc2NhbGVfY29sb3JfbWFudWFsKHZhbHVlcyA9IGNoX2NvbHMsIGd1aWRlPSJub25lIikgK3hsYWIoIkNocm9tb3NvbWUiKSsKICAgIHlsYWIoImxvZzEwKFAtdmFsdWUpIikrCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gbG9nMTAoMC4wNSksIGxpbmV0eXBlID0gJ2Rhc2hlZCcsIGNvbG9yID0gJ3JlZCcsIHNpemU9MC4zKSArCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gLWxvZzEwKDAuMDUpLCBsaW5ldHlwZSA9ICdkYXNoZWQnLCBjb2xvciA9ICdyZWQnLCBzaXplPTAuMykrCiAgICAgIGdndGl0bGUoIlAtdmFsdWVzIGZvciBjaGFuZ2VzIGluIFRoZXRhIikrCiAgICBzY2FsZV94X2NvbnRpbnVvdXMoYnJlYWtzPWNocm1heCRtaWRkbGUsIGxhYmVscz0xOjI2KSsKICAgIHRoZW1lX2J3KCkKZ2dzYXZlKCcuLi9PdXRwdXQvUGkvU2h1ZmZsZS9DaGFuZ2VzX2luX3RoZXRhVy5zaXRlc2h1ZmZsZS5wLXZhbHVlc19QV1MucG5nJywgd2lkdGggPSA3LjUsIGhlaWdodCA9MTEsIHVuaXRzID0gJ2luJywgZHBpID0gMzAwKQoKI3Bsb3QgVGFqYW1hJ3MgRCBwLXZhbHVlIHZzLiBwb3NpdGlvbgpnZ3Bsb3QoRGF0YWxsLCBhZXMoUE9TZ2VuLCAtbG9nMTAodERkLnApKnNpZ24odERkKSwgY29sb3IgPSBDaHJvbW8pKSArIAogICAgZ2VvbV9wb2ludChzaXplID0gMC40LCBhbHBoYSA9IDAuMykgKwogICAgZmFjZXRfd3JhcCh+cG9wLCBuY29sID0gMSkgKwogICAgc2NhbGVfY29sb3JfbWFudWFsKHZhbHVlcyA9IGNoX2NvbHMsIGd1aWRlPSJub25lIikgK3hsYWIoIkNocm9tb3NvbWUiKSsKICAgIHlsYWIoImxvZzEwKFAtdmFsdWUpIikrCiAgICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSBsb2cxMCgwLjA1KSwgbGluZXR5cGUgPSAnZGFzaGVkJywgIGNvbG9yID0gJ3JlZCcsIHNpemU9MC4zKSArCiAgICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSAtbG9nMTAoMC4wNSksIGxpbmV0eXBlID0gJ2Rhc2hlZCcsIGNvbG9yID0gJ3JlZCcsIHNpemU9MC4zKSsKICAgIGdndGl0bGUoIlAtdmFsdWVzIGZvciBjaGFuZ2VzIGluIFRhamltYSdzIEQiKSsKICAgIHNjYWxlX3hfY29udGludW91cyhicmVha3M9Y2hybWF4JG1pZGRsZSwgbGFiZWxzPTE6MjYpKwogICAgdGhlbWVfYncoKQpnZ3NhdmUoJy4uL091dHB1dC9QaS9TaHVmZmxlL0NoYW5nZXNfaW5fVGFqaW1hRC5zaXRlc2h1ZmZsZS5wLXZhbHVlc19QV1MucG5nJywgd2lkdGggPSA3LjUsIGhlaWdodCA9MTEsIHVuaXRzID0gJ2luJywgZHBpID0gMzAwKQoKYGBgCgohW10oLi4vT3V0cHV0L1BpL1NodWZmbGUvQ2hhbmdlc19pbl9QaS5zaXRlc2h1ZmZsZS5wLXZhbHVlc19QV1MucG5nKXt3aWR0aD03MCV9CgoKIVtdKC4uL091dHB1dC9QaS9TaHVmZmxlL0NoYW5nZXNfaW5fdGhldGFXLnNpdGVzaHVmZmxlLnAtdmFsdWVzX1BXUy5wbmcpe3dpZHRoPTcwJX0KCiFbXSguLi9PdXRwdXQvUGkvU2h1ZmZsZS9DaGFuZ2VzX2luX1RhamltYUQuc2l0ZXNodWZmbGUucC12YWx1ZXNfUFdTLnBuZyl7d2lkdGg9NzAlfQoKIyMgRmluZCB0aGUgcmVnaW9ucyB3aXRoIHNpZ25maWljYW50IFAtdmFsdWVzIChPdXRsaWVyIHJlZ2lvbnMpICAKCiMjIyBQaQpgYGB7ciBldmFsPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQojIE91dGxpZXJzCgpQaV9vdXRsaWVyczwtRGF0YWxsW3RQZC5wIDwgMC4wNSxdClRoZXRhX291dGxpZXJzPC1EYXRhbGxbdFdkLnAgPCAwLjA1LF0KVGFqaW1hRF9vdXRsaWVyczwtRGF0YWxsW3REZC5wIDwgMC4wNSxdICNubyBvdXRsaWVycwoKIyBTdW1tYXJ5IHBlciBjaHJvbW9zb21lIHBlciBwb3B1bGF0aW9uCnBpPC1kYXRhLmZyYW1lKHRhYmxlKFBpX291dGxpZXJzJHBvcCwgUGlfb3V0bGllcnMkQ2hyb21vKSkKdGhlPC1kYXRhLmZyYW1lKHRhYmxlKFRoZXRhX291dGxpZXJzJHBvcCwgVGhldGFfb3V0bGllcnMkQ2hyb21vKSkKRDwtZGF0YS5mcmFtZSh0YWJsZShUYWppbWFEX291dGxpZXJzJHBvcCwgVGFqaW1hRF9vdXRsaWVycyRDaHJvbW8pKQoKIyBTdW1tYXJ5IHBlciBwb3B1bGF0aW9uCnBpMjwtZGF0YS5mcmFtZSh0YWJsZShQaV9vdXRsaWVycyRwb3ApKQp0aGUyPC1kYXRhLmZyYW1lKHRhYmxlKFRoZXRhX291dGxpZXJzJHBvcCkpCkRzMjwtZGF0YS5mcmFtZSh0YWJsZShUYWppbWFEX291dGxpZXJzJHBvcCkpCgoKI3Bsb3QgUFdTOTEtOTYsIDk2LTA3LCBhbmQgMDctMTcKeXJzPC1jKCJQV1M5MV9QV1M5NiIsIlBXUzk2X1BXUzA3IiwiUFdTMDdfUFdTMTciKQpjb2wzPC1icmV3ZXIucGFsKDUsIlB1UmQiKVtjKDIsMyw1KV0KZGl2MTwtZGl2ZXJnaW5nX2hjbCg2LCBwYWxldHRlPSJCbHVlLVJlZCIpCiNyZXZlcnNlIHRoZSBjb2xvcgpkaXYyPC1yZXYoZGl2MSkKCnBpczwtcGlbcGkkVmFyMSAlaW4lIHlycyxdCnBpcyRWYXIxPC1mYWN0b3IocGlzJFZhcjEsIGxldmVscz15cnMpCmdncGxvdChwaXMsIGFlcyh4PVZhcjIsIHk9RnJlcSwgZmlsbD1WYXIxKSkrCiAgICBnZW9tX2JhcihzdGF0PSJpZGVudGl0eSIscG9zaXRpb249cG9zaXRpb25fZG9kZ2Uod2lkdGg9MC44KSkrCiAgICBzY2FsZV9maWxsX21hbnVhbCh2YWx1ZXM9Y29sMykrCiAgICB0aGVtZV9jbGFzc2ljKCkrCiAgICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZT00NSwgaGp1c3Q9MSksIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfYmxhbmsoKSkrCiAgICBnZ3RpdGxlKGV4cHJlc3Npb24ocGFzdGUoIkNoYW5nZXMgaW4gIiwgcGkpKSkrCiAgICB4bGFiKCcnKSt5bGFiKCdOdW1iZXIgb2YgcmVnaW9ucyB3aXRoIFA8MC4wNScpCmdnc2F2ZSgiLi4vT3V0cHV0L1BpL1NodWZmbGUvUGlfc2lnbmlmaWNhbnRfcGVyQ2hyb21fcGVyUG9wLnBuZyIsIHdpZHRoID0gOCwgaGVpZ2h0ID0gNCwgZHBpPTMwMCkgCgojIFBlciBwb3B1bGF0aW9uIGNvbXBhcmlzb24KZ2dwbG90KHBpMiwgYWVzKHg9VmFyMSwgeT1GcmVxLCBmaWxsPVZhcjEpKSsKICAgIGdlb21fYmFyKHN0YXQ9ImlkZW50aXR5Iixwb3NpdGlvbj1wb3NpdGlvbl9kb2RnZSh3aWR0aD0wLjgpKSsKICAgIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcz1kaXYxKSsKICAgIHRoZW1lX2NsYXNzaWMoKSsKICAgIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlPTQ1LCBoanVzdD0xKSwgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKSsKICAgIGdndGl0bGUoZXhwcmVzc2lvbihwYXN0ZSgiQ2hhbmdlcyBpbiAiLCBwaSkpKSsKICAgIHhsYWIoJycpK3lsYWIoJ051bWJlciBvZiByZWdpb25zIHdpdGggUDwwLjA1JykKZ2dzYXZlKCIuLi9PdXRwdXQvUGkvU2h1ZmZsZS9QaV9zaWduaWZpY2FudF9wZXJDb21wYXJpc29uLnBuZyIsIHdpZHRoID0gNSwgaGVpZ2h0ID0gNCwgZHBpPTMwMCkgCgojb3JkZXJlZApwaTM8LXBpMltvcmRlcihwaTIkRnJlcSwgZGVjcmVhc2luZyA9IFQpLF0KcGkzJFZhcjE8LWZhY3RvcihwaTMkVmFyMSwgbGV2ZWxzPXBhc3RlMCh1bmlxdWUocGkzJFZhcjEpKSkKZ2dwbG90KHBpMywgYWVzKHg9VmFyMSwgeT1GcmVxLCBmaWxsPVZhcjEpKSsKICAgIGdlb21fYmFyKHN0YXQ9ImlkZW50aXR5Iixwb3NpdGlvbj1wb3NpdGlvbl9kb2RnZSh3aWR0aD0wLjgpKSsKICAgIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcz1kaXYyKSsKICAgIHRoZW1lX2NsYXNzaWMoKSsKICAgIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlPTQ1LCBoanVzdD0xKSwgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKSsKICAgIGdndGl0bGUoZXhwcmVzc2lvbihwYXN0ZSgiQ2hhbmdlcyBpbiAiLCBwaSwgIiAoc2h1ZmZsZSkiKSkpKwogICAgeGxhYignJykreWxhYignTnVtYmVyIG9mIHJlZ2lvbnMgd2l0aCBQPDAuMDUnKQpnZ3NhdmUoIi4uL091dHB1dC9QaS9TaHVmZmxlL1BpX3NpZ25pZmljYW50X3BlckNvbXBhcmlzb25fb3JkZXJlZC5wbmciLCB3aWR0aCA9IDUsIGhlaWdodCA9IDQsIGRwaT0zMDApIAoKI0Fzc2lnbiBjb2xvcnMgdG8gdGhlIGNvbXBhcmlzb24gZ3JvdXAKbmFtZXMoZGl2Mik8LSB1bmlxdWUocGkzJFZhcjEpCgoKYGBgCiFbXSguLi9PdXRwdXQvUGkvU2h1ZmZsZS9QaV9zaWduaWZpY2FudF9wZXJDb21wYXJpc29uX29yZGVyZWQucG5nKSAgCgoKIVtdKC4uL091dHB1dC9QaS9TaHVmZmxlL1BpX3NpZ25pZmljYW50X3BlckNocm9tX3BlclBvcC5wbmcpICAKCgojIyMgVGhldGEKCmBgYHtyIGV2YWw9RkFMU0UsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CgojVGhldGVhCnRoczwtdGhlW3RoZSRWYXIxICVpbiUgeXJzLF0KdGhzJFZhcjE8LWZhY3Rvcih0aHMkVmFyMSwgbGV2ZWxzPXlycykKZ2dwbG90KHRocywgYWVzKHg9VmFyMiwgeT1GcmVxLCBmaWxsPVZhcjEpKSsKICAgIGdlb21fYmFyKHN0YXQ9ImlkZW50aXR5Iixwb3NpdGlvbj1wb3NpdGlvbl9kb2RnZSh3aWR0aD0wLjgpKSsKICAgIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcz1jb2wzKSsKICAgIHRoZW1lX2NsYXNzaWMoKSsKICAgIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlPTQ1LCBoanVzdD0xKSwgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKSsKICAgIGdndGl0bGUocGFzdGUwKCJDaGFuZ2VzIGluIHRoZXRhIikpKwogICAgeGxhYignJykreWxhYignTnVtYmVyIG9mIHJlZ2lvbnMgd2l0aCBQPDAuMDUnKQpnZ3NhdmUoIi4uL091dHB1dC9QaS9TaHVmZmxlL1RoZXRhX3NpZ25pZmljYW50X3BlckNocm9tX3BlclBvcC5wbmciLCB3aWR0aCA9IDUsIGhlaWdodCA9IDMsIGRwaT0zMDApIAoKZ2dwbG90KHRoZTIsIGFlcyh4PVZhcjEsIHk9RnJlcSwgZmlsbD1WYXIxKSkrCiAgICBnZW9tX2JhcihzdGF0PSJpZGVudGl0eSIscG9zaXRpb249cG9zaXRpb25fZG9kZ2Uod2lkdGg9MC44KSkrCiAgICBzY2FsZV9maWxsX21hbnVhbCh2YWx1ZXM9ZGl2MSkrCiAgICB0aGVtZV9jbGFzc2ljKCkrCiAgICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZT00NSwgaGp1c3Q9MSksIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfYmxhbmsoKSkrCiAgICBnZ3RpdGxlKHBhc3RlKCJDaGFuZ2VzIGluIHRoZXRhIikpKwogICAgeGxhYignJykreWxhYignTnVtYmVyIG9mIHJlZ2lvbnMgd2l0aCBQPDAuMDUnKQpnZ3NhdmUoIi4uL091dHB1dC9QaS9TaHVmZmxlL1RoZXRhX3NpZ25pZmljYW50X3BlckNvbXBhcmlzb24ucG5nIiwgd2lkdGggPSA1LCBoZWlnaHQgPSA0LCBkcGk9MzAwKSAKCgp0aGUzPC10aGUyW29yZGVyKHRoZTIkRnJlcSwgZGVjcmVhc2luZyA9IFQpLF0KdGhlMyRWYXIxPC1mYWN0b3IodGhlMyRWYXIxLCBsZXZlbHM9cGFzdGUwKHVuaXF1ZSh0aGUzJFZhcjEpKSkKCmdncGxvdCh0aGUzLCBhZXMoeD1WYXIxLCB5PUZyZXEsIGZpbGw9VmFyMSkpKwogICAgZ2VvbV9iYXIoc3RhdD0iaWRlbnRpdHkiLHBvc2l0aW9uPXBvc2l0aW9uX2RvZGdlKHdpZHRoPTAuOCkpKwogICAgc2NhbGVfZmlsbF9tYW51YWwobmFtZSA9ICJWYXIxIix2YWx1ZXMgPSBkaXYyKSsKICAgIHRoZW1lX2NsYXNzaWMoKSsKICAgIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlPTQ1LCBoanVzdD0xKSwgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKSsKICAgIGdndGl0bGUocGFzdGUoIkNoYW5nZXMgaW4gdGhldGEiKSkrCiAgICB4bGFiKCcnKSt5bGFiKCdOdW1iZXIgb2YgcmVnaW9ucyB3aXRoIFA8MC4wNScpCmdnc2F2ZSgiLi4vT3V0cHV0L1BpL1NodWZmbGUvVGhldGFfc2lnbmlmaWNhbnRfcGVyQ29tcGFyaXNvbl9vcmRlcmVkLnBuZyIsIHdpZHRoID0gNSwgaGVpZ2h0ID0gNCwgZHBpPTMwMCkgCgoKCmBgYAohW10oLi4vT3V0cHV0L1BpL1NodWZmbGUvVGhldGFfc2lnbmlmaWNhbnRfcGVyQ29tcGFyaXNvbi5wbmcpCgohW10oLi4vT3V0cHV0L1BpL1NodWZmbGUvVGhldGFfc2lnbmlmaWNhbnRfcGVyQ2hyb21fcGVyUG9wLnBuZykKCgojIyMgVGFqaW1hJ3MgRCAtIE5vIHJlZ2lvbnMgd2l0aCBQIDwgMC4wNSAgCmBgYHtyIGV2YWw9RkFMU0UsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CgojIFRhamltYSdzIEQKZ2dwbG90KERzMiwgYWVzKHg9VmFyMSwgeT1GcmVxLCBmaWxsPVZhcjEpKSsKICAgIGdlb21fYmFyKHN0YXQ9ImlkZW50aXR5Iixwb3NpdGlvbj1wb3NpdGlvbl9kb2RnZSh3aWR0aD0wLjgpKSsKICAgIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcz1kaXYxKSsKICAgIHRoZW1lX2NsYXNzaWMoKSsKICAgIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlPTQ1LCBoanVzdD0xKSwgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKSsKICAgIGdndGl0bGUocGFzdGUwKCJDaGFuZ2VzIGluIFRhamltYSdzIEQiKSkrCiAgICB4bGFiKCcnKSt5bGFiKCdOdW1iZXIgb2YgcmVnaW9ucyB3aXRoIFA8MC4wNScpCmdnc2F2ZSgiLi4vT3V0cHV0L1BpL1NodWZmbGUvVGFqaW1hRF9zaWduaWZpY2FudF9wZXJDb21wYXJpc29uLnBuZyIsIHdpZHRoID0gNSwgaGVpZ2h0ID0gNCwgZHBpPTMwMCkgCgoKRDI8LURbRCRWYXIxICVpbiUgeXJzLF0KRDIkVmFyMTwtZmFjdG9yKEQyJFZhcjEsIGxldmVscz15cnMpCmdncGxvdChEMiwgYWVzKHg9VmFyMiwgeT1GcmVxLCBmaWxsPVZhcjEpKSsKICAgIGdlb21fYmFyKHN0YXQ9ImlkZW50aXR5Iixwb3NpdGlvbj1wb3NpdGlvbl9kb2RnZSh3aWR0aD0wLjgpKSsKICAgIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcz1jb2wzKSsKICAgIHRoZW1lX2NsYXNzaWMoKSsKICAgIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlPTQ1LCBoanVzdD0xKSwgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKSsKICAgICBnZ3RpdGxlKHBhc3RlMCgiQ2hhbmdlcyBpbiBUYWppbWEncyBEIikpKwogICAgeGxhYignJykreWxhYignTnVtYmVyIG9mIHJlZ2lvbnMgd2l0aCBQPjAuMDUnKQojZ2dzYXZlKCIuLi9PdXRwdXQvUGkvU2h1ZmZsZS9UYWppbWFEX3NpZ25pZmljYW50X3BlckNocm9tX3BlclBvcC5wbmciLCB3aWR0aCA9IDgsIGhlaWdodCA9IDQsIGRwaT0zMDApIAoKYGBgCgotIE1vc3QgZGlmZmVyZW5jZXMgZXhpc3QgYmV0d2VlbiAxOTk2IGFuZCAyMDA3Ci0gQ2hyMjUgaGFzIHRoZSBtb3N0IHNpZ25pZmljYW50IHJlZ2lvbnMgZm9yIGNoYW5nZXMgaW4gUGkgYW5kIFRoZXRhCgo8YnI+Cjxicj4KCiMgRnN0IGdlbm9tZS1zY2FuIHVzaW5nIGEgbnVsbCBtb2RlbAoqIGZyb20gUGluc2t5IGV0IGFsLiAyMDIxICBodHRwczovL2dpdGh1Yi5jb20vcGluc2t5bGFiL2NvZEV2b2wKKiBDcmVhdGUgYSBudWxsIGRpc3RyaWJ1dGlvbiBvZiBleHBlY3RlZCBtYXhpbXVtIEZzdAoqIERlZmluZSBhIGdlbm9tZS13aWRlIHAtdmFsdWUgZm9yIGVhY2ggcmVnaW9uIGFzIChyKzEpLyhuKzEpLCB3aGVyZSByIHdhcyB0aGUgbnVtYmVyIG9mIG51bGwgc2ltdWxhdGlvbnMgd2l0aCBtYXhpbXVtIGNoYW5nZSBncmVhdGVyIHRoYW4gb3IgZXF1YWwgdG8gdGhlIGVtcGlyaWNhbCB2YWx1ZSBhbmQgbiB3YXMgdGhlIG51bWJlciBvZiBzaHVmZmxlcwoKIyMgU2h1ZmZsZSBGc3QgYW5kIGNyZWF0ZSBhIG51bGwgbW9kZWwgb2YgbWF4IEZzdAoqIEJhc2VkIG9uIGFuZ3NkIHBlcnNpdGUgRnN0IG91dHB1dAoKMS4gQWZ0ZXIgcnVubmluZyBhbmdzZCAocmVhbFNGUyksIHByaW50IGEgcGVyc2l0ZSBGc3QgZmlsZSAKYGBge2Jhc2h9CiNydW4gRnN0UFdTcHJpbnQuc2gKCi9ob21lL2phbWNnaXJyL2FwcHMvYW5nc2QvbWlzYy9yZWFsU0ZTIGZzdCBwcmludCAgL2hvbWUva3Rpc3QvcGgvZGF0YS9hbmdzZC9hbmFseXNpcy9mc3RfUFdTOTFfUFdTMDdfcGVyc2l0ZV9tYWYwMC5mc3QuaWR4ID4gL2hvbWUva3Rpc3QvcGgvZGF0YS9hbmdzZC9hbmFseXNpcy9mc3RfUFdTOTFfUFdTMDdfcGVyc2l0ZV9tYWYwMC50eHQKYmd6aXAgL2hvbWUva3Rpc3QvcGgvZGF0YS9hbmdzZC9hbmFseXNpcy9mc3RfUFdTOTFfUFdTMDdfcGVyc2l0ZV9tYWYwMC50eHQKCmBgYAoKMi4gUnVuIEZzdF9zaHVmZmxlX3B3cy5SIGF0IEZhcm0gKHNsdXJtIHNjcmlwdDogYW5nc2RfZnN0X3NpdGVzaHVmZmxlX251bGwuc2gpICAKCmBgYHtyIGV2YWw9RkFMU0UsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CiMgRnJvbSBodHRwczovL2dpdGh1Yi5jb20vcGluc2t5bGFiL2NvZEV2b2wgYW5nc2RfZnN0X3NpdGVzaHVmZmxlX251bGwucgoKIyBzaHVmZmxlIEFOR1NEIHBlcnNpdGUgRlNUIEEgYW5kIEIgdmFsdWVzIGFjcm9zcyBzaXRlcyBhbmQgY2FsY3VsYXRlIHdpbmRvd2VkIEZTVCB0byBnZXQgYSBudWxsIGRpc3RyaWJ1dGlvbiBvZiBtYXggZ2Vub21lLXdpZGUgRlNUCiMgbmVlZCBmc3QgcGVyc2l0ZSBvdXRwdXQgZnJvbSBhbmdzZCBmc3RfKl9wZXJzaXRlX21hZjAwLnR4dC5neiBmaWxlcwoKIyB0byBydW4gb24gZmFybQoKCiMjIyMjIFIgY29kZSAgIyMjIyMjIyAKCiMgcGFyYW1ldGVycwp3aW5zeiA8LSA1MDAwMCAjIHdpbmRvdyBzaXplCndpbnN0cCA8LSAxMDAwMCAjIHdpbmRvdyBzdGVwCm5yZXAgPC0gMTAwMCAjIG51bWJlciBvZiByZXNodWZmbGVzCm1pbmxvY2kgPC0gMTAgIyBtaW5pbXVtIG51bWJlciBvZiBsb2NpIHBlciB3aW5kb3cgdG8gY29uc2lkZXIKCgojIGxvYWQgZnVuY3Rpb25zCnJlcXVpcmUoZGF0YS50YWJsZSkKCgojIyMjIyMjIyMjIyMjCiMgUHJlcCBkYXRhCiMjIyMjIyMjIyMjIyMKIyBsb2FkIGZzdCBBL0IgZGF0YQojY2FuIDwtIGZyZWFkKCdhbmFseXNpcy9DYW5fNDAuQ2FuXzE0LmZzdC5BQi5neicpCiNzZXRuYW1lcyhjYW4sIGMoJ0NIUk9NJywgJ1BPUycsICdBJywgJ0InKSkKCnB3czkxOTYgPC0gZnJlYWQoJy9ob21lL2t0aXN0L3BoL2RhdGEvYW5nc2QvYW5hbHlzaXMvZnN0X1BXUzkxX1BXUzk2X3BlcnNpdGVfbWFmMDAudHh0Lmd6JykKc2V0bmFtZXMocHdzOTE5NiwgYygnQ0hST00nLCAnUE9TJywgJ0EnLCAnQicpKQpwd3M5MTA3IDwtIGZyZWFkKCcvaG9tZS9rdGlzdC9waC9kYXRhL2FuZ3NkL2FuYWx5c2lzL2ZzdF9QV1M5MV9QV1MwN19wZXJzaXRlX21hZjAwLnR4dC5neicpCnNldG5hbWVzKHB3czkxMDcsIGMoJ0NIUk9NJywgJ1BPUycsICdBJywgJ0InKSkKcHdzOTExNyA8LSBmcmVhZCgnL2hvbWUva3Rpc3QvcGgvZGF0YS9hbmdzZC9hbmFseXNpcy9mc3RfUFdTOTFfUFdTMTdfcGVyc2l0ZV9tYWYwMC50eHQuZ3onKQpzZXRuYW1lcyhwd3M5MTE3LCBjKCdDSFJPTScsICdQT1MnLCAnQScsICdCJykpCnB3czk2MDcgPC0gZnJlYWQoJy9ob21lL2t0aXN0L3BoL2RhdGEvYW5nc2QvYW5hbHlzaXMvZnN0X1BXUzk2X1BXUzA3X3BlcnNpdGVfbWFmMDAudHh0Lmd6JykKc2V0bmFtZXMocHdzOTYwNywgYygnQ0hST00nLCAnUE9TJywgJ0EnLCAnQicpKQpwd3M5NjE3IDwtIGZyZWFkKCcvaG9tZS9rdGlzdC9waC9kYXRhL2FuZ3NkL2FuYWx5c2lzL2ZzdF9QV1M5Nl9QV1MxN19wZXJzaXRlX21hZjAwLnR4dC5neicpCnNldG5hbWVzKHB3czk2MTcsIGMoJ0NIUk9NJywgJ1BPUycsICdBJywgJ0InKSkKcHdzMDcxNyA8LSBmcmVhZCgnL2hvbWUva3Rpc3QvcGgvZGF0YS9hbmdzZC9hbmFseXNpcy9mc3RfUFdTMDdfUFdTMTdfcGVyc2l0ZV9tYWYwMC50eHQuZ3onKQpzZXRuYW1lcyhwd3MwNzE3LCBjKCdDSFJPTScsICdQT1MnLCAnQScsICdCJykpCgpQV1M8LWxpc3QoKQpQV1NbWzFdXTwtcHdzOTE5NgpQV1NbWzJdXTwtcHdzOTEwNwpQV1NbWzNdXTwtcHdzOTExNwpQV1NbWzRdXTwtcHdzOTYwNwpQV1NbWzVdXTwtcHdzOTYxNyAKUFdTW1s2XV08LXB3czA3MTcgCiAgICAKIAojIGNyZWF0ZSBuZXcgY29sdW1ucyBhcyBpbmRpY2VzIGZvciB3aW5kb3dzCgpmb3IgKGkgaW4gMTogbGVuZ3RoKFBXUykpewogICAgZGY8LVBXU1tbaV1dCiAgICBmb3IoaiBpbiAxOih3aW5zei93aW5zdHApKXsKICAgICAgICBkZlssIChwYXN0ZTAoJ3dpbicsIGopKSA6PSBmbG9vcigoUE9TIC0gKGotMSkqd2luc3RwKS93aW5zeikqd2luc3ogKyB3aW5zei8yICsgKGotMSkqd2luc3RwXQogICAgfQogICAgUFdTW1tpXV08LWRmCn0KCiMgbWFyayB3aW5kb3dzIHdpdGggPCBtaW5sb2NpIGZvciByZW1vdmFsCnJlbSA8LSByZXAoMCwgNikgIyBudW1iZXIgb2Ygd2luZG93cyByZW1vdmVkIGZvciBlYWNoIG9mIHRoZSA2IGNvbXBhcmlzb25zCgpmb3IoaSBpbiAxOiBsZW5ndGgoUFdTKSl7CiAgICBwdzwtUFdTW1tpXV0KICAgIGZvcihqIGluIDE6KHdpbnN6L3dpbnN0cCkpewogICAgICAgIHB3d2luIDwtIHB3WywgLihuc25wcyA9IGxlbmd0aChQT1MpKSwgYnkgPSAuKHdpbiA9IGdldChwYXN0ZTAoJ3dpbicsIGopKSldICMgY2FsYyBudW0gc25wcyBwZXIgd2luZG93CiAgICAgICAgcmVtWzFdIDwtIHJlbVsxXSArIHB3d2luWywgc3VtKG5zbnBzIDwgbWlubG9jaSldICMgcmVjb3JkIG51bWJlciB0byBiZSByZW1vdmVkCiAgICAgICAgcHd3aW5bLCAocGFzdGUwKCd3aW4nLCBqLCAna2VlcCcpKSA6PSAxXSAjIGNyZWF0ZSBjb2wgdG8gbWFyayB3aGljaCB3aW5kb3dzIHRvIGtlZXAKICAgICAgICBwd3dpbltuc25wcyA8IG1pbmxvY2ksIChwYXN0ZTAoJ3dpbicsIGosICdrZWVwJykpIDo9IDBdICMgbWFyayB3aW5kb3dzIHRvIHJlbW92ZQogICAgICAgIHB3d2luWywgbnNucHMgOj0gTlVMTF0gIyBkcm9wIGNvbHVtbgogICAgICAgIHNldG5hbWVzKHB3d2luLCAid2luIiwgcGFzdGUwKCd3aW4nLCBqKSkgIyBjaGFuZ2UgY29sIG5hbWUKICAgICAgICBwdyA8LSBtZXJnZShwdywgcHd3aW4sIGJ5ID0gcGFzdGUwKCd3aW4nLCBqKSwgYWxsLnggPSBUUlVFKSAjIG1lcmdlIGtlZXBlciBjb2wgYmFjayB0byBmdWxsIGRhdGFzZXQKICAgIH0KICAgIFBXU1tbaV1dPC1wdwp9CgpyZW0gIyBudW1iZXIgb2Ygd2luZG93cyByZW1vdmVkIGZvciBlYWNoIGNvbXBhcmlzb24KCgoKCiMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIwojIHNodWZmbGUgYW5kIHJlY2FsYyB3aW5kb3dlZCBGU1QKIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjCmNvbG5tcyA8LSBjKCdDSFJPTScsICdQT1MnLCBwYXN0ZTAoJ3dpbicsIDE6KHdpbnN6L3dpbnN0cCkpLCBwYXN0ZTAoJ3dpbicsIDE6KHdpbnN6L3dpbnN0cCksICdrZWVwJykpICMgbGlzdCBvZiBjb2x1bW4gbmFtZXMgd2Ugd2FudCBvdXQgb2YgdGhlIGJhc2UgZGF0YS50YWJsZQoKIyBQV1MwNy0xNwpwb3BzPC1jKCJQV1M5MS45NiIsIlBXUzkxLjA3IiwiUFdTOTEuMTciLCJQV1M5Ni4wNyIsIlBXUzk2LjE3IiwiUFdTMDcuMTciKQoKCmZvcihwIGluIDE6bGVuZ3RoKFBXUykpewogICAgcHJpbnQocGFzdGUwKCdTdGFydGluZyAnLCBwb3BzW3BdKSkKCiAgICBwdzwtUFdTW1twXV0KICAgIGZvcihpIGluIDE6bnJlcCl7CiAgICAgICAgY2F0KGkpOyBjYXQoJyAnKQogICAgICAgICMgY3JlYXRlIG5ldyBkYXRhc2V0CiAgICAgICAgaW5kcyA8LSBzYW1wbGUoMTpucm93KHB3KSwgbnJvdyhwdyksIHJlcGxhY2UgPSBGQUxTRSkKICAgICAgICB0ZW1wIDwtIGNiaW5kKHB3WywgLi5jb2xubXNdLCBwd1tpbmRzLCAuKEEsIEIpXSkgIyBzaHVmZmxlIEZTVHMgYWNyb3NzIHBvc2l0aW9ucwogICAgICAgIAogICAgICAgICMgY2FsYyBmc3QgZm9yIGVhY2ggd2luZG93IHRvIGtlZXAKICAgICAgICBmb3IoaiBpbiAxOih3aW5zei93aW5zdHApKXsKICAgICAgICAgICAgdGVtcDIgPC0gdGVtcFtnZXQocGFzdGUwKCd3aW4nLCBqLCAna2VlcCcpKSA9PSAxLCBdICMgdHJpbSB0byB3aW5kb3dzIHRvIGtlZXAuIGNhbid0IGNvbWJpbmUgd2l0aCBuZXh0IGxpbmUgZm9yIHNvbWUgcmVhc29uLgogICAgICAgICAgICBpZihqID09MSkgdGVtcGZzdHMgPC0gdGVtcDJbLCAuKGZzdCA9IHN1bShBKS9zdW0oQikpLCBieSA9IC4oQ0hST00sIFBPUyA9IGdldChwYXN0ZTAoJ3dpbicsIGopKSldCiAgICAgICAgICAgIGlmKGogPiAxKSB0ZW1wZnN0cyA8LSByYmluZCh0ZW1wZnN0cywgdGVtcDJbLCAuKGZzdCA9IHN1bShBKS9zdW0oQikpLCBieSA9IC4oQ0hST00sIFBPUyA9IGdldChwYXN0ZTAoJ3dpbicsIGopKSldKQogICAgICAgIH0KICAgICAgICAKICAgICAgICAjIHNhdmUgdGhlIG1heCB3aW5kb3dlZCBmc3QKICAgICAgICAjIGV4Y2x1ZGUgd2luZG93cyB3aXRoIG5lZ2F0aXZlIG1pZHBvaW50cwogICAgICAgIGlmKGkgPT0gMSkgbWF4ZnN0IDwtIHRlbXBmc3RzW1BPUyA+IDAsIG1heChmc3QsIG5hLnJtID0gVFJVRSldCQogICAgICAgIGlmKGkgPiAxKSBtYXhmc3QgPC0gYyhtYXhmc3QsIHRlbXBmc3RzW1BPUyA+IDAsIG1heChmc3QsIG5hLnJtID0gVFJVRSldKQogICAgfQoKICAgIHByaW50KHBhc3RlKCdNYXg6JywgbWF4KG1heGZzdCwgbmEucm0gPSBUUlVFKSwgJzsgOTV0aDonLCBxdWFudGlsZShtYXhmc3QsIHByb2IgPSAwLjk1LCBuYS5ybSA9IFRSVUUpKSkKICAgIAogICAgd3JpdGUuY3N2KG1heGZzdCwgZ3pmaWxlKHBhc3RlMCgnL2hvbWUva3Rpc3QvcGgvZGF0YS9hbmdzZC9hbmFseXNpcy8nLHBvcHNbcF0sICdfc2l0ZXNodXVmZmxlLmNzdi5neicpKSwgcm93Lm5hbWVzID0gRkFMU0UpCiAgICBybShtYXhmc3QpCn0KCmBgYAoKCiMjIEFzc2VzcyB0aGUgcmVzdWx0cyBmcm9tIHNpdGUgcmVzaHVmZmxlZCBGc3QgCk5lZWQgCi0gc2l0ZWh1ZmZsZSByZXN1bHRzICpmc3Rfc2l0ZXNodXVmZmxlLmNzdi5negotIHdpbmRvdy1iYXNlZCBGc3QgcmVzdWx0cyBmcm9tIGFuZ3NkICpfNTBrV2luZG93X21hZjAwCi0gcGVyc2l0ZSBGc3QgZmlsZXMgYWZ0ZXIgTEQgcHJ1bmluZyAod2l0aCBBQiBpbmZvKSAqX3BlcnNpdGVfbWFmMDBfbGR0cmltLmNzdi5neiAKCmBgYHtyIGV2YWw9RkFMU0UsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CiNmcm9tIGh0dHBzOi8vZ2l0aHViLmNvbS9waW5za3lsYWIvY29kRXZvbAoKIyBwYXJhbWV0ZXJzCm1pbmxvY2kgPC0gMTAgCndpbnN6IDwtIDUwMDAwICMgd2luZG93IHNpemUKd2luc3RwIDwtIDEwMDAwICMgd2luZG93IHN0ZXAKCiMgbG9hZCBmdW5jdGlvbnMKcmVxdWlyZShkYXRhLnRhYmxlKQpyZXF1aXJlKGdncGxvdDIpCgpjYWxjcCA8LSBmdW5jdGlvbihmc3QsIG51bGwpIHJldHVybigoc3VtKG51bGwgPiBmc3QpKzEpLyhsZW5ndGgobnVsbCkrMSkpICMgZXF1YXRpb24gZnJvbSBOb3J0aCBldCBhbC4gMjAwMiBBbSBKIEh1bSBHZW4KCiMgcmVhZCBpbiBhbmQgcHJlcCBkYXRhCnllYXJzPC1jKCJQV1M5MSIsIlBXUzk2IiwiUFdTMDciLCJQV1MxNyIpCmNvbWI8LXQoY29tYm4oeWVhcnMsIDIpKQoKIyBjb250aW51b3VzIG51Y2xlb3RpZGUgcG9zaXRpb24gZm9yIHRoZSB3aG9sZSBnZW5vbWUKY2hybWF4IDwtIGZyZWFkKCcuLi9EYXRhL25ld192Y2YvY2hyX3NpemVzLmJlZCcpCmNocm1heDwtY2hybWF4WywtMl0KY29sbmFtZXMoY2hybWF4KTwtYygiY2hyIiwgImxlbiIpCmNocm1heCRzdGFydDwtYygwLGN1bXN1bShjaHJtYXgkbGVuKVsxOihucm93KGNocm1heCktMSldKQpzZXRrZXkoY2hybWF4LCBjaHIpCgpwcnVuZWRGc3RfZ3c8LWRhdGEuZnJhbWUocG9wMT1jb21iWywxXSxwb3AyPWNvbWJbLDJdKQojZm9yIChwIGluIDE6IG5yb3coY29tYikpewogICAgZm9yIChwIGluIDE6NSl7CiAgICAjZ2Vub21lLXdpZGUgbWF4IEZzdCBmcm9tIHJlc2h1ZmZsaW5nICh1bmxpbmtlZCBzaXRlcykKICAgIG51bGw8LWZyZWFkKHBhc3RlMCgnRGF0YS9zaHVmZmxlLycsIGNvbWJbcCwxXSwiLiIsY29tYltwLDJdLCdfZnN0X3NpdGVzaHV1ZmZsZS5jc3YuZ3onKSkKICAgICN3aW5kb3ctYmFzZWQgRnN0IGZyb20gYW5nc2QgICAgICAgICAgICAKICAgIGZzdHdpbiA8LSBmcmVhZChwYXN0ZTAoJ0RhdGEvbmV3X3ZjZi9hbmdzZC9mcm9tVkNGLzJEL2ZzdF8nLGNvbWJbcCwxXSwiXyIsY29tYltwLDJdLCdfNTBrV2luZG93X21hZjAwJyksIGhlYWRlciA9IEZBTFNFLCBjb2wubmFtZXMgPSBjKCdyZWdpb24nLCAnY2hyJywgJ21pZFBvcycsICdOc2l0ZXMnLCAnZnN0JykpICMgYWxsIHNpdGVzCiAgICAjcGVyc2l0ZSBmc3QgKEFCKSAtcHJ1bmVkIHNpdGVzCiAgICBBQiA8LSBmcmVhZChwYXN0ZTAoJ0RhdGEvc2h1ZmZsZS9mc3RfJyxjb21iW3AsMV0sIl8iLGNvbWJbcCwyXSwnX3BlcnNpdGVfbWFmMDBfbGR0cmltLmNzdi5neicpKQogICAgCiAgICAjIG1lcmdlIG51Y2xlb3RpZGUgcG9zaXRpb24gaW50byB0aGUgZnJlcXVlbmN5IGZpbGVzCiAgICBzZXRrZXkoQUIsIENIUk9NKQogICAgQUIgPC0gQUJbY2hybWF4WywgLihDSFJPTSA9IGNociwgc3RhcnQpXSwgXQogICAgQUJbLCBwb3NnZW4gOj0gUE9TICsgc3RhcnRdCiAgICBBQlssc3RhcnQgOj0gTlVMTF0KICAgIAogICAgIyBDYWxjIGdlbm9tZS13aWRlIEZTVAogICAgcHJ1bmVkRnN0X2d3JGd3RnN0W3BdPC1BQlshaXMubmEoQSksIHN1bShBKS9zdW0oQildCiAgICAKICAgICMgQ2FsYyB3aW5kb3dlZCBGU1QKICAgICMgY3JlYXRlIG5ldyBjb2x1bW5zIGFzIGluZGljZXMgZm9yIHdpbmRvd3MKICAgIGZvcihqIGluIDE6KHdpbnN6L3dpbnN0cCkpewogICAgICAgIEFCWywgKHBhc3RlMCgnd2luJywgaikpIDo9IGZsb29yKChQT1MgLSAoai0xKSp3aW5zdHApL3dpbnN6KSp3aW5zeiArIHdpbnN6LzIgKyAoai0xKSp3aW5zdHBdCiAgICB9CiAgICAKICAgICMgY2FsYyBmc3QgYW5kICMgc25wcyBwZXIgd2luZG93CiAgICBmb3IoaiBpbiAxOih3aW5zei93aW5zdHApKXsKICAgICAgICBpZihqID09MSl7CiAgICAgICAgICAgIGZzdHdpbiA8LSBBQlssIC4oZnN0ID0gc3VtKEEpL3N1bShCKSwgbmxvY2kgPSBsZW5ndGgoUE9TKSksIGJ5ID0gLihDSFJPTSwgbWlkUG9zID0gZ2V0KHBhc3RlMCgnd2luJywgaikpKV0KICAgICAgICB9IAogICAgICAgIGlmKGogPiAxKXsKICAgICAgICAgICAgZnN0d2luIDwtIHJiaW5kKGZzdHdpbiwgQUJbLCAuKGZzdCA9IHN1bShBKS9zdW0oQiksIG5sb2NpID0gbGVuZ3RoKFBPUykpLCBieSA9IC4oQ0hST00sIG1pZFBvcyA9IGdldChwYXN0ZTAoJ3dpbicsIGopKSldKQogICAgICAgIH0gCiAgICB9CiAgICAKICAgICMjIG51bGwgbW9kZWwgc3RhdHMKICAgIGZzdGF0czwtbnVsbFssIC4obWF4ID0gbWF4KHgpLCB1OTUgPSBxdWFudGlsZSh4LCBwcm9icyA9IDAuOTUpKV0KICAgIHBydW5lZEZzdF9ndyRudWxsLkZzdG1heFtwXTwtZnN0YXRzJG1heFsxXQogICAgcHJ1bmVkRnN0X2d3JG51bGwuRnN0Ljk1cVtwXTwtZnN0YXRzJHU5NVsxXQogICAgCiAgICAjQ2FsY3VsYXRlIHAtdmFsdWVzCiAgICBwdmFsPC1mc3R3aW5bLCBwIDo9IGNhbGNwKGZzdCwgbnVsbCR4KSwgYnkgPSAuKENIUk9NLCBtaWRQb3MpXQogICAgCiAgICAjY29tYmluZWQgdGhlIGRhdGFzZXRzCiAgICBwYWlyPC1wYXN0ZTAoY29tYltwLDFdLCJfIixjb21iW3AsMl0pCiAgICBmc3R3aW5bLCBwb3AgOj0gcGFpciBdCiAgICB3cml0ZS5jc3YoZnN0d2luLCBwYXN0ZTAoIi4uL091dHB1dC9Gc3QvZnN0X3NpdGVzaHVmZmxlLnBhaXJ3aXNlXyIsIHBhaXIsIi5jc3YiKSkKfQoKd3JpdGUuY3N2KHBydW5lZEZzdF9ndywgIi4uL091dHB1dC9Gc3QvZnN0X3NpdGVzaHVmZmxlX3B3c19zdW1tYXJ5LmNzdiIpCgpkYXRhPC1kYXRhLmZyYW1lKCkKZm9yKGkgaW4gMTogbnJvdyhjb21iKSl7CiAgICBwYWlyPC1wYXN0ZTAoY29tYltpLDFdLCJfIixjb21iW2ksMl0pCiAgICBkZjwtZnJlYWQocGFzdGUwKCIuLi9PdXRwdXQvRnN0L2ZzdF9zaXRlc2h1ZmZsZS5wYWlyd2lzZV8iLCBwYWlyLCIuY3N2IikpCiAgICBkZjwtZGZbIWlzLm5hKGZzdCkgJiBtaWRQb3MgPiAwLCBdCiAgICBkZlssY2g6PSBhcy5pbnRlZ2VyKGdzdWIoImNociIsJycsIENIUk9NKSldCiAgICBzZXRrZXkoZGYsQ0hST00pCiAgICBkZjwtZGZbY2hybWF4WywgLihDSFJPTSA9IGNociwgc3RhcnQpXSwgXQogICAgZGZbLHBvcy5ndzo9IG1pZFBvcytzdGFydF0KICAgICNzZXRrZXkoZGYsIGNoLCBtaWRQb3MpCiAgICBkYXRhPC1yYmluZChkYXRhLCBkZikKfQoKd3JpdGUuY3N2KGRhdGEsIGZpbGUgPSBnemZpbGUoJy4uL091dHB1dC9Gc3QvZnN0X3NpdGVzaHVmZmxlLmFuZ3NkLlBXUy5jc3YuZ3onKSwgcm93Lm5hbWVzID0gRkFMU0UpCmBgYAoKCiMjIFBsb3QgdGhlIHJlc3VsdHMgICAgICAKYGBge3IgZXZhbD1GQUxTRSwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KCmRhdGE8LWZyZWFkKCcuLi9PdXRwdXQvRnN0L2ZzdF9zaXRlc2h1ZmZsZS5hbmdzZC5QV1MuY3N2Lmd6JykKCiNwb3AgYXMgZmFjdG9yIApkYXRhWywgcG9wIDo9IGZhY3Rvcihwb3AsIGxldmVscyA9IGMoIlBXUzkxX1BXUzk2IiwiUFdTOTFfUFdTMDciLCJQV1M5MV9QV1MxNyIsIlBXUzk2X1BXUzA3IiwiUFdTOTZfUFdTMTciLCJQV1MwN19QV1MxNyIpKV0KCiMgcGxvdCBwLXZhbHVlIHZzLiBwb3NpdGlvbgp0d29fY29sczwtcmVwKGMoInN0ZWVsYmx1ZSIsImxpZ2h0Ymx1ZSIpLCB0aW1lcz0xMykKbWlubG9jaSA8LSAxMCAjIG1pbmltdW0gbnVtYmVyIG9mIGxvY2kgcGVyIHdpbmRvdyB0byBjb25zaWRlcgoKZ2dwbG90KGRhdGFbbmxvY2kgPj0gbWlubG9jaSwgXSwgYWVzKHBvcy5ndywgLWxvZzEwKHApLCBjb2xvciA9IGZhY3RvcihjaCkpKSArIAogICAgZ2VvbV9wb2ludChzaXplID0gMC4yLCBhbHBoYSA9IDAuNSkgKwogICAgZmFjZXRfd3JhcCh+cG9wLCBuY29sID0gMSkgKwogICAgc2NhbGVfY29sb3JfbWFudWFsKHZhbHVlcyA9IHR3b19jb2xzLCBndWlkZT0nbm9uZScpICsKICAgIGdlb21faGxpbmUoeWludGVyY2VwdCA9IC1sb2cxMCgwLjA1KSwgbGluZXR5cGUgPSAnZGFzaGVkJywgY29sb3IgPSAnZ3JleScpKwogICAgdGhlbWVfYncoKSsKICAgIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF9ibGFuaygpKSsKICAgIHhsYWIoIkdlbm9tZSIpCmdnc2F2ZSgiLi4vT3V0cHV0L0ZzdC9mc3Quc2l0ZXNodWZmbGUucF92c19wb3MuUFdTLnBuZyIsIHdpZHRoID0gNy41LCBoZWlnaHQgPSA3LjUsZHBpID0gMzAwKQpgYGAKIVtdKC4uL091dHB1dC9Gc3QvZnN0LnNpdGVzaHVmZmxlLnBfdnNfcG9zLlBXUy5wbmcpCgotIE51bWJlciBvZiBsb2NpIHBlciB3aW5kb3cgdnMuIEZzdCBwLXZhbHVlcwoKYGBge3IgZXZhbD1GQUxTRSwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KIyBwbG90IHAtdmFsdWUgdnMuIG5sb2NpIHRvIGNoZWNrIGlmIGFueSBiaWFzZXMgZXhpc3QKZ2dwbG90KGRhdGEsIGFlcyhubG9jaSwgLWxvZzEwKHApLCBjb2xvciA9IHBvcCkpICsgCiAgICBnZW9tX3BvaW50KHNpemUgPSAwLjYsIGFscGhhID0gMC41KSArCiAgICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSAtbG9nMTAoMC4wNSksIGxpbmV0eXBlID0gJ2Rhc2hlZCcsIGNvbG9yID0gJ2dyZXknKSArIAogICAgc2NhbGVfeF9sb2cxMCgpK3RoZW1lX2J3KCkrdGhlbWUobGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKSsKICAgIGd1aWRlcyhjb2xvciA9IGd1aWRlX2xlZ2VuZChvdmVycmlkZS5hZXMgPSBsaXN0KHNpemUgPSAzLCBhbHBoYT0xKSkpCmdnc2F2ZSgiLi4vT3V0cHV0L0ZzdC9mc3Quc2l0ZXNodWZmbGUucF92c19ubG9jaS5QV1MucG5nIiwgd2lkdGggPSA4LCBoZWlnaHQgPTQsZHBpID0gMzAwKQpgYGAKCiFbXSguLi9PdXRwdXQvRnN0L2ZzdC5zaXRlc2h1ZmZsZS5wX3ZzX25sb2NpLlBXUy5wbmcpe3dpZHRoPTcwJX0KCiMjIFNpZ25pZmljYW50IG91dGxpZXIgcmVnaW9ucyAgCgojIyMgcC12YWx1ZSA8IDAuMDUKCmBgYHtyIGVjaG89VFJVRSwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0Kb3V0bGllcnM8LWRhdGFbcCA8IDAuMDUsIF0Kd3JpdGUuY3N2KG91dGxpZXJzLCAiLi4vT3V0cHV0L0ZzdC9Gc3RfbnVsbHNodWZmbGluZ19vdXRsaWVyc19QV1MuY3N2IikKCmNvdW50czwtZGF0YS50YWJsZSh0YWJsZShvdXRsaWVycyRwb3ApKQpnZ3Bsb3QoY291bnRzLCBhZXMoeD1WMSwgeT1OKSkrCiAgICBnZW9tX2JhcihzdGF0PSJpZGVudGl0eSIsIGZpbGw9Y29sc1s0XSkrCiAgICB4bGFiKCcnKSt5bGFiKCcnKSsKICAgIHRoZW1lX2NsYXNzaWMoKSsKICAgIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlPTQ1LCBoanVzdD0xKSkKZ2dzYXZlKCIuLi9PdXRwdXQvRnN0L0ZzdF9vdXRsaWVyc19ieUNvbXBhcmlzb24ucG5nIiwgd2lkdGggPSA0LCBoZWlnaHQgPSAzLCBkcGk9MzAwICkKCgpjb3VudHMkVjE8LWZhY3Rvcihjb3VudHMkVjEsIGxldmVscz1jKCJQV1M5Nl9QV1MwNyIsIlBXUzA3X1BXUzE3IiwiUFdTOTFfUFdTOTYiLCAiUFdTOTFfUFdTMDciLCAiUFdTOTFfUFdTMTciLCJQV1M5Nl9QV1MxNyIpKQpnZ3Bsb3QoY291bnRzLCBhZXMoeD1WMSwgeT1OLCBmaWxsPVYxKSkrCiAgICBnZW9tX2JhcihzdGF0PSJpZGVudGl0eSIpKwogICAgc2NhbGVfZmlsbF9tYW51YWwobmFtZSA9ICJWIix2YWx1ZXMgPSBkaXYyKSsKICAgIHhsYWIoJycpK3lsYWIoJ051bWJlciBvZiByZWdpb25zIHdpdGggUDwwLjA1JykrCiAgICB0aGVtZV9jbGFzc2ljKCkrZ2d0aXRsZSgiRnN0IChzaHVmZmxlKSIpKwogICAgdGhlbWUoYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGU9NDUsIGhqdXN0PTEpLCBsZWdlbmQudGl0bGUgPSBlbGVtZW50X2JsYW5rKCkpCmdnc2F2ZSgiLi4vT3V0cHV0L0ZzdC9Gc3Rfb3V0bGllcnNfYnlDb21wYXJpc29uX29yZGVyZWQucG5nIiwgd2lkdGggPSA0LCBoZWlnaHQgPSAzLCBkcGk9MzAwICkKCgoKI291dGxpZXJzPC1yZWFkLmNzdigiLi4vT3V0cHV0L0ZzdC9Gc3RfbnVsbHNodWZmbGluZ19vdXRsaWVyc19QV1MuY3N2Iiwgcm93Lm5hbWVzID0gMSkKa25pdHI6OmthYmxlKG91dGxpZXJzWyxjKDIsMyw0LDUsNiw3LDgpXSwgY2FwdGlvbj0iT3V0bGllciByZWdpb25zIikKCiNPdXRsaWVyIHJlZ2lvbnMKI0NIUk9NCW1pZFBvcwlmc3QJbmxvY2kJcAlwb3AJY2gKI2NocjI1CTgyNTAwMAkwLjEyMTA5OTcJMTM4CTAuMDAzOTk2CVBXUzk2X1BXUzA3CTI1CiNjaHIyNQk4MzUwMDAJMC4wODE4ODAyCTI3OQkwLjAyMDk3OQlQV1M5Nl9QV1MwNwkyNQojY2hyMjUJODQ1MDAwCTAuMDcyMTkyNgkzMTEJMC4wMzM5NjYJUFdTOTZfUFdTMDcJMjUKI2NocjI1CTgwNTAwMAkwLjEyNDI2NTIJMTY2CTAuMDAyOTk3CVBXUzk2X1BXUzA3CTI1CiNjaHIyNQk4MTUwMDAJMC4xMzc3MjE4CTE1NQkwLjAwMjk5NwlQV1M5Nl9QV1MwNwkyNQojY2hyMTkJMTM3NTAwMAkwLjEyNTE2MDIJNDIJMC4wMDM5OTYJUFdTMDdfUFdTMTcJMTkKCmBgYAohW10oLi4vT3V0cHV0L0ZzdC9Gc3Rfb3V0bGllcnNfYnlDb21wYXJpc29uLnBuZykgCgoKIyMjIFRvcCAxJSBvdXRsaWVycyAmIFA8IDAuMgoqIHJlbGF4aW5nIHRoZSBjdXRvZmYKYGBge3J9CiNkYXRhPC1mcmVhZCgnLi4vT3V0cHV0L0ZzdC9mc3Rfc2l0ZXNodWZmbGUuYW5nc2QuUFdTLmNzdi5neicpCmRhdDwtZGF0YVtvcmRlcihwKSxdCnRvcDE8LWRhdFsxOihucm93KGRhdCkqMC4wMSksXQp0b3AxPC10b3AxW3A8MC4yLF0KCmNvdW50czI8LWRhdGEudGFibGUodGFibGUodG9wMSRwb3ApKQpnZ3Bsb3QoY291bnRzMiwgYWVzKHg9VjEsIHk9TikpKwogICAgZ2VvbV9iYXIoc3RhdD0iaWRlbnRpdHkiLCBmaWxsPWNvbHNbNF0pKwogICAgeGxhYignJykreWxhYignJykrCiAgICB0aGVtZV9jbGFzc2ljKCkrCiAgICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZT00NSwgaGp1c3Q9MSkpCmdnc2F2ZSgiLi4vT3V0cHV0L0ZzdC9Gc3Rfb3V0bGllcnNfYnlDb21wYXJpc29uX3RvcDFwZXJjZW50LnBuZyIsIHdpZHRoID0gNCwgaGVpZ2h0ID0gMywgZHBpPTMwMCApCgoKCgprbml0cjo6a2FibGUodG9wMVssYygyLDMsNCw1LDYsNyw4KV0pCmBgYAohW10oLi4vT3V0cHV0L0ZzdC9Gc3Rfb3V0bGllcnNfYnlDb21wYXJpc29uX3RvcDFwZXJjZW50LnBuZykgIAoKIVtdKC4uL091dHB1dC9Gc3QvcGNhbmdzZF9zY2FuX2hpZ2hQc2l0ZXMucG5nKSAgCiogKipOTyBvdmVybGFwcGluZyByZWdpb25zIGZyb20gUENBbmdzZCBzZWxlY3Rpb24gc2NhbioqCgoKCgoKCmBgYHtyIGV2YWw9RkFMU0UsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CgpgYGAKCmBgYHtyIGV2YWw9RkFMU0UsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CgpgYGAK